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ABSTRACT 


\ 

This  report  covers  work  performed  under  Contract 
AF23( 601) -4009  with  the  Aeronautical  Chart  and  Information 
Center,  St.  Louis,  Missour^M^v  the  Hawaii  Institute  of  Geophysics 

of  the  University  of  Hawaii.  I ^presents  the  results  of  theor- 

r  c\  '.'t'J 

etical  studies  fracasied  oat  ia~oa»der  to  determine  the  combined 
geodetic  and  geophysical  significance  of  current  gravity  re¬ 
duction  techniques.  Special  emphasis  was  given  to  the  problem 
of  interpolation  of  deflections  of  the  vertical  by  gravimetric 
means . 

An  examination  of  the  mathematical  theory  of  physical 
geodesy  revealed  that  the  accuracy  of  the  mathematical  equations 
relating  deflections  of  the  vertical  and  gravity  anomalies  could 
not  be  improved  by  choice  of  gravity  anomaly.  The  choice  of  a 
gravity  anomaly  for  use  in  the  deflection  interpolation  equations 
was  found  to  depend  upon  the  accuracy  with  which  the  anomaly 
could  be  interpolated  between  observation  points  and  the 
relative  amount  of  labor  required  to  carry  out  numerical 
computations . 

The  complete  Bouguer  anomaly  with  geologic  corrections  was 
found  to  be  the  best  anomaly  to  use  for  interpolation  of  gravity 
with  a  reversion  to  the  normal  complete  Bouguer  anomaly  without 
geologic  corrections  for  use  in  the  mathematical  equations. 


TABLE  OF  CONTENTS 


» 


Page 


SECTION  1. 

Introduction 

1 

SECTION  2. 

Basic  Potential  Theory 

5 

The  Boundary  Value  Problem 

5 

Development  in  Spherical  Harmonics 

7 

SECTION  3- 

Classical  Geodetic  Theory 

13 

Deflections  of  the  Vertical  in  Classical  Theory 

26 

SECTION  ij-. 

The  New  Geodetic  Theory 

30 

SECTION  5. 

Theoretical  Aspects  of  Anomaly  Selection 

^7 

Classical  Theory 

bl 

The  New  Theory 

55 

SECTION  6. 

Problems  of  Interpolation  of  Anomalies 

60 

The  Bouguer  Anomaly 

63 

Isostatic  Anomalies 

70 

Model  Earth  Anomalies 

73 

Rudzki  Anomalies 

75 

Summary 

77 

SECTION  7. 

Formulation  for  Solution 

8o 

BIBLIOGRAPHY 

90 

l 


SECTION  1 


INTRODUCTION 

This  report  covers  the  results  of  the  Investigation  and  Selection 
Phases  of  the  Alpine  Gravity  Reduction  Project.  As  stated  in  the  technical 
specifications,  the  work  to  be  carried  out  in  the  Investigation  and  Selection 

Phases  was  to:  l)  " - evaluate  current  gravity  reduction  techniques  for  their 

combined  geodetic  and  geophysical  significance”  and  2)  " - select  that  gravity 

reduction  method  which  will  most  accurately  reproduce  the  external  gravity  field.” 
The  over-all  aim  of  this  project  was  to  find  the  best  method  of  utilizing  geologic 
and  geophysical  information  in  determining  gravity  anomalies  to  be  used  in  the 
interpolation  of  deflections  of  the  vertical  in  areas  of  rugged  topography. 

In  order  to  determine  what  type  of  gravity  anomaly  would  be  the  best 
to  use  for  interpolation  of  deflections  of  the  vertical  in  practical  problems, 
two  distinct  investigations  were  necessary.  First,  the  mathematical  theory  lead¬ 
ing  to  the  equations  relating  gravity  anomalies  and  deflections  of  the  vertical 
was  examined.  In  this  way  it  was  possible  to  determine  if  the  mathematical  theory 
itself  led  to  any  special  preference  for  a  particular  type  of  gravity  anomaly  and 
if  controls  were  placed  upon  the  type  of  anomaly  which  could  be  used  by  the  nature 
of  the  mathematical  formulation.  Second,  since  gravity  is  actually  known  only  at 
a  limited  number  of  points,  the  various  types  of  gravity  anomalies  were  studied 
to  determine  their  relative  value  in  terms  of  accuracy  of  interpolation  between 
points  of  observation. 

Since  the  equations  relating  gravity  anomalies  and  deflections  of  the 
vertical  are  obtained  by  solving  a  boundary  value  problem,  a  study  of  the  mathe¬ 
matical  theory  is  a  study  of  this  boundary  value  problem.  The  boundary  value 
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problem  can  be  stated  in  words  as  follows:  "Given  the  gravitational  potential 
and  the  magnitude  of  the  force  of  gravity  everywhere  on  the  bounding  surface  of 
a  gravitating  body,  determine  the  shape  of  the  bounding  surface  and  the  direction 
of  the  force  of  gravity  on  it."  For  a  body  such  as  the  earth,  whose  bounding 
surface  is  highly  irregular,  some  approximations  must  be  made  in  order  to  make 
the  problem  tractable.  There  have  been  two  modes  of  approach  to  simplification 
of  the  problem.  In  the  sections  which  follow,  these  have  been  labeled  "The 
Classical  Geodetic  Theory"  and  "The  New  Geodetic  Theory".  These  labels  were  not 
chosen  to  have  any  connotation  as  to  accuracy  but  were  chosen  simply  because  the 
approach  labeled  "The  Classical  Geodetic  Theory"  was  the  original  approach  used 
in  solving  the  problem  while  the  approach  labeled  "The  New  Geodetic  Theory"  has 
been  developed  only  recently. 

The  "Classical  Theory"  is  developed  in  Section  3.  Section  4  gives  the 
development  of  the  "New  Theory".  Section  2  contains  background  material  which 
gives  the  development  of  certain  relations  used  in  Sections  3  and  4  and  is  primarily 
to  be  used  for  reference  when  reading  these  sections. 

In  the  "classical"  approach,  as  developed  in  Section  3,  a  solution  is 
obtained  for  the  boundary  value  problem  under  the  simplifying  assumption  that 
the  bounding  surface  is  an  equipotential  surface.  In  this  case,  the  type  of 
gravity  anomaly  used  is  controlled  by  the  necessity  for  computing  the  effect  of 
the  theoretical  transfer  of  all  masses  of  the  actual  earth  which  lie  outside  the 
geoid  to  some  point  within  the  geoid.  In  Section  5,  the  various  types  of  anomalies 
which  will  accomplish  such  a  mass  transfer  are  examined. 

In  the  "new"  approach,  as  developed  in  Section  4,  the  free-air  anomaly 
arises  naturally  in  the  development.  However,  as  shown  in  Section  5,  this  natural 
appearance  of  the  free-air  anomaly  is  simply  a  result  of  the  theoretical  model 
chosen.  With  other  choices  of  theoretical  models,  other  types  of  gravity  anomalies 
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can  be  used.  The  conclusions  reached  in  Sections  3>  4  and  5  are:  l)  the  "New 
Theory"  is  capable  of  extensions  to  greater  accuracy  than  is  the  "Classical  Theory". 
The  increased  accuracy  is  of  limited  importance  for  computation  of  height  dif¬ 
ferences  but  becomes  significant  for  deflection  computations  in  rugged  terrain. 

2)  In  using  the  "Classical  Theory"  it  is  difficult  to  determine  which  type  of 
anomaly  gives  the  most  accurate  result  since  the  answer  appears  to  depend  in  a 
complex  way  upon  the  difference  between  the  actual  density  distributions  above 
sea  level  and  the  actual  gravity  gradients  and  the  values  assumed  for  these 
quantities  when  carrying  out  the  calculations.  3)  In  using  the  "New  Theory",  it  is 
found  that  gravity  anomalies  computed  using  any  number  of  density  models  will  give 
equally  accurate  results  so  long  as  all  of  the  mass  of  each  model  is  assumed  to 
lie  within  the  earth’s  surface. 

Since,  in  the  "New  Theory"  the  accuracy  of  the  results  is  theoretically 
independent  of  the  type  of  gravity  anomaly  used,  the  choice  of  anomaly  depends,  in 
the  case  of  practical  computations,  upon:  l)  the  relative  accuracy  with  which  the 
various  types  of  anomalies  can  be  interpolated,  and  2)  the  amount  of  work  required 
to  compute  a  particular  type  of  anomaly  and  carry  out  computations  using  it.  In 
Section  6  the  various  types  of  anomalies  are  examined  to  determine  which  can  best 
satisfy  these  two  criterias.  It  is  here  in  connection  with  the  accuracy  of  inter¬ 
polation  that  the  question  of  the  degree  to  which  the  model  densities  represent 
the  actual  densities  of  the  earth  becomes  important  when  using  the  "New  Theory". 

From  the  point  of  view  of  interpolation  accuracy,  one  would  like  to  use  a  density 
model  which  approaches  as  nearly  as  possible  to  the  actual  density  distribution 
of  the  earth  above  some  fixed  depth,  this  depth  being  dependent  upon  the  distance 
between  points  at  which  gravity  has  been  observed.  In  the  case  of  deflection 
interpolation  computations  which  would  not  normally  be  attempted  without  a  reason¬ 
ably  close -spaced  net  of  observations,  this  leads  to  some  type  of  Bouguer  anomaly 


-4- 


vith  geologic  corrections  applied  down  to  whatever  depth  is  necessary  to  cause 
the  resultant  anomaly  to  vary  smoothly  between  observation  points. 

Using  such  geologically  corrected  gravity  anomalies  directly  in  the 
deflection  formulae  was  found  to  lead  to  complicated  computation  procedures. 

The  compromise  finally  adopted  was,  therefore,  to  use  the  geologic  control  for 
gravity  interpolation  but  to  revert  to  normal  complete  Bouguer  anomalies  for 
computation . 

The  method  of  solution  chosen  was  one  developed  by  Pellinen  (1962), 
which  utilizes  terrain  corrected  free-air  anomalies.  The  theory  leading  to  the 
final  deflection  formula  using  this  method  is  given  in  Section  7.  The  theory  is 
presented  in  a  form  which  allows  two  contour  maps — one  of  complete  Bouguer  anom¬ 
aly  and  one  of  elevation — to  be  used  rather  them  a  single  map  of  terrain  corrected 
free-air  anomalies.  This  method  of  procedure  has  the  advantage  that  the  anomaly 
used,  the  Complete  Bouguer  Anomaly,  is  a  smoothly  varying  function  while  the 
terrain  corrected  free-air  anomaly  is  nearly  linearly  related  to  local  elevation 
changes  and  therefore  changes  in  a  complex  manner  in  regions  of  rugged  terrain. 

Detailed  computational  procedures  using  the  theory  of  Section  7  are 
given  in  Part  II  of  the  final  report  on  this  contract  in  conjunction  with  the  re¬ 
sults  obtained  from  applying  this  theory  to  actual  problems  in  the  Test  and 
Application  Phases  of  the  contract. 
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SECTION  2 

Basic  Potential  Theory 

The  purpose  of  this  section  Is  to  suunarlze  those  developments  In 
potential  theory  which  are  utilized  in  the  development  of  the  various  theories  of 
physical  geodesy.  Only  the  basic  outlines  will  be  given  here  with  the  aim  of 
providing  an  easily  accessable  source  of  reference  for  some  of  the  basic  formulae. 
For  a  more  detailed  and  mathematically  rigorous  development  one  may  refer  to 
MacMillan  (1938).  Hie  results  here  are  largely  a  summary  of  the  results  presented 
there. 


2.1  The  Boundary  Value  Problem 

We  shall  begin  with  the  well  known  Green* s  equation  which  states:  If  9 
and  ♦  are  any  two  functions  of  the  Cartesian  coordinates  x,y,  and  z,  and  if  their 
first  derivatives  are  continuous  and  the  functions  possess  second  derivatives 
within  a  volume  V  and  on  its  surface  S,  then  the  following  integral  relation  holds 
between  9  and  ♦  . 

2 

9  -  9  a  fjdv  -  j.  if  «  -  9  ^)dS 

(2-1) 


Jy  (f  A2  «p  .  «p  A2  t)dV  -/,<♦§!-*  ^)dS 


2  a2  .a2  a2 

where  A  =  the  Laplacian  operator  ~2  ***  — J  **■  — 2 

BX  BY  BZ 

-5.  =  the  derivative  in  the  direction  of  the  outward  normal 
dn  to  the  surface  S,  and 

and  J*B  =  integrals  over  the  volume  and  surface  respectively. 

For  the  development  of  this  theory  the  reader  is  referred  to  MacMillan  sections  53* 
54,  and  55. 

For  use  In  potential  theory.  Green's  formula  is  utilized  as  follows.  If 
the  functions 9  and  f  of  (2-l)  satisfy  Laplace’s  equation  A2  9  .  o  and  A2  f  .  0 
within  the  volume  V  bounded  by  the  surface  S  then  (2-l)  becomes 


(2-2) 
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It  is  noted  that  if  S.TLC  are  the  coordinates  of  a  variable  point  in  V  and  if 


x,y,  and  z  are  the  coordinates  of  any  point  P  then  the  function 


t.  [<« . 


X)2  +  (T|  -  y)2 


+  (C  -  *)2]' 


(2-3) 


satisfies  the  conditions  necessary  for  use  as  one  of  the  functions  in  (2-2) 
everywhere  in  the  volume  V  provided  P  lies  outside  the  volume.  If  P  lay  within 
the  volume  V  itself  problems  would  arise  where  the  variable  point  approached  the 
point  P.  Thus  for  P  outside  V  we  can  use  ♦  =  ~  in  (2-2)  and  have  the  following 


equation 


?  II  •  *  si  (|>  ds  ■  0 


(2-4) 


For  our  purposes  we  will  let  the  volume  V  of  equation  (2-l)  be  all  of  the 
space  lying  outside  the  gravitating  body  which  we  are  studying  except  for  a  small 
sphere  around  the  point  P  which  is  assumed  to  lie  outside  the  gravitating  body. 

Then  the  surface  of  the  volume  V  is  the  surface  of  the  body  S’  and  the  surface  of 
the  small  sphere  around  P  which  we  shall  call  S**.  We  pause  here  to  re-emphasize 
that  the  volume  V  is  all  the  area  outside  the  gravitating  body  and  the  small  sphere 


enclosing  P.  Then  since  P  lies  outside  V,  we  can  write  (2-4)  as 

J s ’  [l  II  -  <p  s!  (?}]  ds'  +  -Ts"  Il  -  *  5!  (?>j 


dS"  -  0 


However,  it  can  be  shown  (See  MacMillan,  1958)  that 


j's"  [?  In  -  *  5!  (?)  ds"  ■ 


(2-5) 


(2-6) 


where  the  subscript  P  indicates  the  function  cp  is  to  be  evaluated  at  the 


point  P. 


Then  (2-5)  becomes 


Jr  [1 

[p  si  - 


^  <5~n 


(2-7) 
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Stated  in  words  (2-7)  says  that  if  a  function  ♦  satisfies  the  conditions: 

1.  9  and  its  first  derivatives  are  single  valued  and  continuous  within 

the  volume  V  and  on  its  boundary  S. 

2.  the  second  derivatives  of  9  exist 

3.  9  satisfies  the  equation  of  Laplace 

within  a  volume  which  is  bound  by  a  closed  surface  S*  then  the  value  of  9  at  any 
interior  point  of  the  volume  can  be  computed  from  the  values  of  9  and  its  first 
spatial  derivatives  on  S'. 

In  geodetic  theory,  the  function  9  will  be  the  potential  of  some  mass 
distribution  lying  entirely  outside  the  volume  V  (i.e.,  within  the  gravitating 

body  bounded  by  S*)  or  on  the  surface  S*.  Such  a  potential  function  may  be  used 

2 

since  it  satisfies  A  9  *  0  outside  the  gravitating  body.  The  potential  related 
to  centrifugal  force  does  not  satisfy  the  Laplace  equation  and  thus  is  not  a 
function  which  will  satisfy  (2-7) •  It  is  for  this  reason  that  in  the  development 
of  the  theory  of  Section  3  the  potential  T  must  be  considered  as  purely  a  gravita¬ 
tion  potential.  In  this  case  T  satisfies  (2-7)  and  we  can  write 

TP '  I .  •  [?  H  -  T  s!  ds’  (2-8) 

2.2  Development  in  Spherical  Harmonics 

Since  development  in  spherical  harmonics  is  fundamental  in  physical 
geodesy,  two  ways  of  arriving  at  a  spherical  harmonic  representation  of  gravitation¬ 
al  potential  are  outlined  below.  It  is  believed  that  taken  together,  the  two  give 
a  clear  understanding  of  what  spherical  harmonics  imply. 

Let  us  begin  the  first  development  by  considering  the  gravitational 
potential, at  a  point  P  lying  outside  a  given  mass  distribution,  due  to  a  small 
incremental  volume  V,  whose  mass  1b  dM  and  which  is  located  at  a  point  Q.  Using 
the  symbol  U  to  represent  potential  we  have 

kdM  kodV 

dUP  ’  ““  *  (2-9) 
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where:  k  =  gravitational  constant 

dV  ss  incremental  volume  element 
a  =  density 

r  =  distance  from  P  to  Q 

Thus  the  gravitational  potential  at  the  point  P  due  to  the  entire  mass  distribu¬ 
tion  is  given  by 


I?  *k  f  <7(r,,6,,X')  j  t  jflf  .  \» 

up  k  Jv  7TTr;vt7TT7'p7*7TJ  dr  d9  dX 

where  the  coordinates  have  the  meaning  illustrated  in  Figure  (2-l) 
From  Figure  (2-l)  we  see  by  the  law  of  cosines  that 


(2-10) 


r  -  (p"  +  r'2  .  2p r 1  cos  *)1/2  .  p(i  +  (£l) 2  _  2 


1/2 


COS  t  ) 


1/2 


(2-11) 


Thus 


i  -  i  (1  +  (I-!)2  .  2  (f)  coa  ♦ )-  1/2 


(2-12) 


The  binomial  series  expansion  is  given  by 

(1  +  x)n 


,  ,  __  .  n(n-l)  2  n  (n-  1 )  (n-  2 )  3 

l  +  nx  T - j-f —  x  t - if *t - x 


TT 


+  -  -  -  - 


(2-13) 


provided  x<l  .  Thus,  if  in  equation  (2-12)  p>r*  we  can  use  (2-13)  to  expand 

r'  2  r  * 

the  quantity  in  parenthesis  in  (2-12)  with  n  =  -  l/2  and  x  -I  (p-)  -  2  (— )  cos  f] 

r» 

If  we  can  carry  out  such  an  expansion  and  collect  terms  with  common  powers  of  (-^) 
we  get  i  *  -i  |l  +  (--)  coa  t  +  (-^-)2  co®2  ^  -  1/2)  +  0--)3 


(5/2  cos  t  -  3/2  coa  t) + 


} 


(2-14) 


It  can  be  shewn  that  the  function  of  cos  t  appearing  In  the  nth  term  of  the  ex¬ 
pansion  can  be  written 


_ 1 dn(co»  f2  -  l)n 

2nn.'  d  cos 


(2-15) 
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The  functions  generated  by  this  formula  are  called  Legendre  polynomials  or  zonal 
harmonics  and  are  represented  by  the  symbol  Pn  (cos  ♦  ).  Using  the  Pn  symbology 


we  can  therefore  write 


?  “Jo  ^n+I  Pn(C0S  ♦) 


(2-16) 


Thus  taking  note  of  (2-l4)  and  (2-l6)  we  can  write  (2-10)  as 

K  “.So  <“•  ♦>«- i  J.  ”  "■ 

+  iL  J  a  r  1  2  (3/2  cos  *  -  l/2)dV+— (2-17) 

P  V 

It  is  instructive  to  evaluate  the  first  few  integrals  of  the  right  hand 


side  of  equation  (2-17). 


The  first  integral  is  simply 


Jv  adV  =  Mg  _  total  mass  of  the  earth 

The  second  integral  can  be  written 

J*  °  r*  cos  tdV  -  cos  9  J  o  z’dV  +  sin  0 


(2-18) 


cos  X  Jv  tfx'dV  +  si"  0  «in  X  Jv  ay  '  dV  (2-19) 
But  the  coordinates  (x,  y,  z)  of  the  center  of  mass  of  a  body  are  by 


definition 


Jv  CTx'dV  _  Jv  ay'dv  _  J'v°E’dV 
- FT 5  y - TT i  * - m: 


(2-20) 


where  Mg  is  the  mass  of  the  body.  Thus  if  we  choose  the  origin  of  our  coordinate 
system  to  be  the  center  of  mass  of  the  body,  then  x  =  y  =  z  =  0  which  leads  to  the 


result 


J  ox’dv  “  Jv  tfy'dV  *  Jv  a*'dV  ■  0 


(2-21) 


The  third  integral  on  the  right  hand  side  of  (2-17)  can  he  vrittep 
Jv  o"r ' 2  (3/2  cos2  *  -  1/2)  -  (3/2  cos2  8  .  1/2)  Jy  <?  «’2 

-  1/2  (x 1 2  +  y,2)jdv  +  3/2  sin  28  cos  X  Jvx's'dV  +  3/2  sin  28  sin  X 
r  ay'z'dv  +  3/2  sin2  8  sin  2X  J*  Ox'y'dy+  3/2  sin2  8  cos  2X 


J val/2(x'2  -  y'2)dV  (2-22) 
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If  the  z*  axis  of  the  coordinate  system  coincides  with  one  of  the 
principal  axes  of  inertia  of  the  gravitating  body 


J  ax's 'dV  *  J  ay 'x'dv  -  0 


(2-23) 


In  the  case  of  the  earth  we  chose  the  z*  axis  to  coincide  with  the 
rotational  axis  of  the  earth,  which  we  know  is  to  a  very  close  approximation  a 
principal  axis  of  inertia,  so  that  we  can  assume  (2-23)  holds. 

Thus,  taking  all  of  the  above  results  into  account  we  can  write  (2-17)  as 


ft  *  p  [< 


(3/2  cos  8  _  1/2) 


[«* 


2  -  1/  2  (x  '  2  +  y 


dV  + 


3/2  sin20  sin  2X  f  ox'y'dV  +  3/2  sip2  9  cos  2X  Jy  ol/2(x  -  y  )<*Vj 

+  L  ff  J3  ~STI  Pn  <c°8  ♦>«v}  (2'2U) 

The  above  development  shows  that  Mtne  constants  in  spherical  harmonic 
development  of  the  gravitational  field  are  various  moments  and  demonstrates  the 
reason  for  certain  terms  being  zero.  However,  the  above  method  does  not  demon¬ 
strate  how  to  get  the  general  expansion  in  terms  of  9  and  X  .  Briefly  this  is 
arrived  at  by  seeking  a  general  solution  to  Laplaces  equation  in  spherical 
coordinates  (P  ,  9  ,  X  )  i  .  e  . 

i  !  a  ,p2aD.  i  a  .  .  a  auv  .  1  a2w\  0 

p  is?  (“5ir)  +  Tnn  X*  (8ln  0  tb>  +  ’  0  (2-25) 


The  general  solution  is  found  to  be 


0  -  Jo  An  p"  Y  +  K  B  l2_ 

n“°  n  n  n=0  n  n+T 

P 


(2-26) 


where  the  and  Bn  are  constants  and 


Yn  =i=0  “'nl  Fn  (cos  9)  cos  1  X  +ilobni  Pn  (co8  9)  8ln  1  X 

where  ani  and  b^  are  also  constants  and 

P°  (cos  8)  -  __L - 2 — ^  (cos2  8  -  l)n 

n  2%:  dcose" 

PX  (cos  0)  =  (cos2  9-l)i/2  — - — (cos2  9-l)n  , 

n  dcos8n+1  1  *  0 


(2-27) 


(2-28) 
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Since  the  potential  function  must  vanish  at  infinity,  must  be  zero 
and  it  is  conventional  procedure  to  set  all  Bn  equal  to  the  gravitational 
constant  k,then  (2-2 6)  becomes 


U 


Y 

k  2  n 

n*0  ~n+I 


(2-29) 


Equation  (2-29)  is  a  general  equation  expressing  the  potential  of  a  mass 
distribution  in  spherical  harmonics.  In  the  practical  case  in  physical  geodesy 
we  wish  to  express  the  attraction  potentials  of  the  actual  earth  and  of  the 
ellipsoidal  normal  field  in  spherical  harmonics  and  then  to  take  the  difference 
to  get  a  spherical  harmonic  representation  of  the  anomalous  potential  T.  In 
normal  practice  the  center  of  gravity  of  the  actual  earth  is  chosen  as  the  center 
of  the  coordinate  system  to  be  used  and  the  center  of  gravity  of  the  normal  field 
is  chosen  to  coincide  with  the  c.g.  of  the  actual  earth.  In  this  case,  as  has 
been  shown,  the  Y^  term  will  be  zero  for  both  potentials  and  thus  for  their  dif¬ 
ference  T.  It  is  also  normal  to  choose  the  mass  of  the  normal  model  to  equal  the 
mass  of  the  actual  earth.  In  this  case,  the  Y0  terms  will  be  identical  and  will 

cancel  when  computing  T.  Thus  it  is  permissable  to  write  T  as 

00  Y 

T  =  If  T  n 

“-2  ~n+l  (2-30) 

The  symbol  Y^  is  used  to  represent  the  spherical  harmonic  of  degree  n 
when  referring  to  a  number  of  different  bodies  throughout  this  report.  To  avoid 
confusion,  one  should  point  out  that  the  symbol  does  not  refer  to  the  same 
quantity  for  all  bodies.  For  any  particular  body  the  constants  of  (2-27)  have 
particular  values  characteristic  of  that  body  and  related  to  the  moments  of 


the  body. 
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SECTION  3 

Classical  Geodetic  Theory 

The  present  section  discusses  the  determination  of  the  shape  of  the 
bounding  equipotential  surface  of  body  B  according  to  the  theory  usually  em¬ 
ployed  in  classical  geodesy. 

Let  V  be  a  reference  gravitational  plus  rotational  field  whose 
corresponding  potential  field  is  represented  by  U.  In  practice,  for  computation 
of  gravity  anomalies,  the  reference  field  is  taken  to  be  that  given  by  the 
International  Gravity  formula  with  the  gravity  gradient  taken  to  be  a  constant 
.09^06  mgls/ft.  A  number  of  papers  have  been  written  on  more  exact  expressions 
for  the  normal  field.  These  have  been  examined  by  Daugherty  and  Define  (1963) 
who  found  that  for  geodetic  purposes  all  are  sufficiently  accurate. 

If  the  gravitational  plus  rotational  potential  of  the  body  B  is 
represented  by  the  symbol  W,  the  potential  at  any  point  is  given  by 

W  =  U  +  T  (3-1) 

where  T  is  by  definition  the  difference  between  W  and  U. 

Consider  the  figure  below: 
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The  equipotential  surface  W  =  Wp  is  taken  to  be  the  bounding 

equipotential  surface,  S,  of  the  body  B.  P  is  a  point  on  this  surface  S.  Up  is 

the  equipotential  of  the  normal  field  passing  through  the  point.  U  is  the  equi¬ 
ty 

potential  surface  of  the  reference  field  such  that 

Uq  =  Wp  (3-2) 

The  point  Q  is  the  point  on  the  equipotential  surface  Uq  intersected  by  the 
normal  to  the  potential  field  U  which  passes  through  P.  N  is  the  distance  PQ. 

Two  points  should  be  mentioned  here.  First,  there  are  other  methodB 
of  establishing  the  point  Q  on  Uq  to  be  related  to  the  point  P  on  Wp.  These 
lead  to  slightly  different  forms  of  development  leading  to  equation  (3-5)  below, 
but  the  final  results  Eire  essentially  the  same  (See,  for  example,  Jung,  1956). 
Second,  the  use  of  (3-2)  implies  that  the  value  of  Wp  is  known  so  that  one  can 
establish  which  equipotential  surface  of  the  normal  field  Uq  is  equal  to  it. 

In  fact,  for  problems  concerned  with  the  earth,  Wp  is  not  known.  For  a  dis¬ 
cussion  of  the  errors  which  result,  see  Molodenski  et  al  (i960).  We  shall  pro¬ 
ceed  under  the  assumption  that  Wp  is  known. 

From  definition  (3-l)  we  have 

Wp  =  Up  +  Tp  (3-3) 

Now  assume  that  if  the  potential  function  U  is  expanded  about  the  point 
Q  in  a  Taylor  series,  the  value  of  U  at  P  can  be  accurately  enough  represented  by 
the  first  two  terms  of  the  expansion.  Then  we  have 


UP  ■  UQ  -  0Q»  =  -  V 


(3-4) 


where  the  negative  sign  results  from  the  fact  that  the  direction  of  y  is  toward 
Q  from  P.  The  accuracy  of  the  assumption  that  the  first  two  terms  are  sufficient 
depends  upon  the  magnitude  of  the  distance  N.  If  we  consider  that  geoidal  heights 
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seldom  exceed  100  feet,  we  can  get  am  idea  of  the  magnitude  of  the  error 

N  =  100  feet  N2  - 

y  m  980,000  mg Is .  dY 


Then  YqN  =  98  X  106  to  108 


Tv 


10,000 

M  .1  Bgls . / ft . 


m2  «  .lx  10,000  -  1000  -  103 

Thus  for  reasonable  N  values,  the  neglected  terms  are  negligible  (i.e.,  of  the 
order  10"5  of  the  second  term) . 

Then  substituting  (3-4)  into  (3-3)  gives 

Wp 

or  making  use  of  (3-2) 


UQ  -  YQK  +  TP 


K  -  IE 

Q 


(3-5) 


In  classical  geodetic  theory  the  bounding  equipotential  Wp  is  the  co- 
geoid  differing  from  the  geoid  or  sea  level  equipotential  of  the  earth  by  an 
amount  which  depends  upon  the  method  used  to  reduce  the  surface  anomalies  of  the 
earth  to  geoid  level  anomalies  and  the  equipotential  surface  Uq  is  the  inter¬ 
national  ellipsoid.  Then  N  is  the  distance  between  ellipsoid  and  co-geoid.  The 
implications  of  the  reduction  methods  used  in  classical  geodesy  will  be  discussed 
in  Section  5*  For  the  present  assume  gravity  is  known  on  Wp. 

Since  Yq  is  known, only  Tp  must  be  determined  in  order  to  compute  N.  To 
do  this  we  will  establish  a  differential  equation  relating  gravity  on  the  equi¬ 
potential  surface  W  =  Wp  to  Tp  and  then  solve  this  equation. 

Referring  to  (3-l)  we  see  that  it  is  possible  to  write  the  identity 


,5W,  ,3  IK  .  /3T. 

Wp  =  Wp  +  Wp 


(3-6) 
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vhere  the  subscripted  P  indicates  the  derivatives  are  to  be  evaluated  at  the 


point  P. 


(Sn'p  8p 


(3-7) 


the  "observed”  gravity  at  P  and 


(%^)  *  %v>  cos  e  -  Yp  cos  € 


(3-8) 


where  t  is  the  angle  between  the  normals  at  P  to  the  equipotentials  Wp  and  Up, 
i.e.,  it  is  the  "deflection  of  the  vertical"  at  P.  Since  e  seldom  exceeds  1’, 
the  approximation  cos  «  =  1  is  permissable.  Using  this  approximation  and 
substituting  (3-7)  and  (3-8)  into  (3-6)  gives 


gp  ■  YP  +  <«>, 


(3-9) 


Expanding  the  function  y  about  the  point  Q  in  a  Taylor  series  and 
keeping  only  two  terms  in  the  expansion  yields 


VQ  "  (Iv)qN 


(3-io) 


Again  note  that  the  error  here  is  of  the  same  order  as  discussed  after  equation 

(3-k). 

Substituting  for  N  from  (3-5)  gives 


Tp 

WQ  7^ 


Substituting  this  into  (3-9)  gives 


c  =  Y  \  Tp  /BT» 

8P  YQ  -  o v ' q  +  <*n>p 

,  V\  /8Tv  Tp  ,dY* 

(gp  -  YQ>  *  (5-)p  -  y-  (?_)Q 


(3-11) 
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This  formula  is  called  the  boundary  value  problem  of  classical  geodesy. 

It  describes  the  relation  between  gravity  anomalies  and  T  at  points  on  the  boundary 
of  the  body.  For  the  case  of  the  earth  both  U  and  W  are  potentials  of  attraction 
plus  a  rotational  potential.  It  is  usually  assumed  that  the  two  rotational  poten¬ 
tials  are  equal. 

Thus,  T  =  W  -  U  is  a  pure  gravitational  potential  and  is  an  harmonic 
function  outside  the  body  B.  The  assumption  that  the  rotational  potential  of  W 
and  U  are  equal  may  not  be  strictly  accurate.  The  problems  connected  with  this 
possibility  are  discussed  for  example  by  Molodenski,  et  al  (i960).  Normally, 
this  problem  is  ignored. 

If  the  reference  gravitational  field  were  spherical  rather  than  the 

Y  .  kM 

ellipsoidal  normal  field  usually  employed  we  would  have  Q  ~7  and 

*Q 


<H> 


Wq 


(3-13) 


where  k  =  universal  gravitational  constant 

M  =  mass  of  model  usually  assumed  to  be  equal  to  mass  of  earth 
Rq=  distance  from  center  of  spherical  field  to  Q 

Substituting  equations  (3-13)  into  equation  (3-12)  gives 


^8P  ”  R  TP 


Now  if  it  is  also  sufficiently  accurate  to  use  the  approximation 


(3-14) 


/dr »  /3  x  v 

Wp  -  -  w, 


(3-15) 


we  get 


dT, 


(8p  -  V  •  ■  R  Tp  •  (n’p 


(3-16) 
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In  the  derivation  of  Stokes  and  Venlng  Melnesz  equations  the  boundary 
condition  (3-16)  is  used.  Note  that  in  deriving  (3-16)  the  approximations 
(3-13)  and  (3-15)  are  applied  to  the  right  hand  side  of  equation  (3-12)  but  the 
values  of  used  on  the  left  hand  side  of  equation  (3-16)  in  deriving  the  gravity 
anomalies  remain  the  International  Formula  values,  i.e.,  the  ellipsoidal  approxi¬ 
mation  is  retained  on  the  left  hand  side  of  the  equation. 

The  error  involved  in  using  equation  (3-16)  rather  than  (3-12)  would  be 
expected  to  be  of  the  order  of  the  flattening  of  the  earth.  The  result  of 
proceeding  from  (3-12)  retaining  the  assumption  of  an  ellipsoidal  normal  field  and 
obtaining  a  result  correct  to  the  order  of  the  square  of  the  flattening  has  been 
studied  most  intensively  by  Zagrebin  (1956)  and  by  Bjerhammar  (1962)  who  cor¬ 
rected  some  errors  in  Zagrebin' s  final  formulae.  The  results  for  computing  N  are 
very  complicated  and  have  received  little  usage.  Molodenski  et  al  (i960)  outline 
earlier  work  of  Molodenski  and  arrive  at  a  somewhat  simpler  formulation  which  he 
extends  to  deflection  of  the  vertical  computations  in  the  classical  theory.  We 
shall  not  at  present  pursue  this  refinement  further  insofar  as  it  should  not  have 
any  significant  effect  on  deflection  interpolation  problems. 

The  next  step  in  the  normal  procedure  is  to  determine  a  formula  for  T 
which  will  be  correct  everywhere  outside  of  and  on  the  surface,  S,  of  the  body  B. 
The  value  of  T  anywhere  on  or  outside  a  gravitating  body  can  be  expressed  by  the 
use  of  Green's  formula  (See  Section  2).  Letting  P  be  any  point  on  or  outside  the 
surface  S,  we  have 

Tp  -  l/4n  [i  -  *£5  <|>lds 

(3-17) 

where  r  =  distance  from  P  to  variable  point  on  surface  S  and  n  is  the  normal  to 


the  surface  S  directed  inward. 
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Cons  ider  the  situation  below: 


r'  =  distance  between  0  and  incremental  surface  element  dS 
r  =  distance  from  P  to  incremental  surface  element  dS 
p  =  distance  from  P  to  0 
t  =  angle  between  r*  and  P 

Assume  that  the  point  P  lies  sufficiently  far  from  the  origin  that, 
for  all  elements  dS  of  the  surface,  P  >  r  *.  Then,  as  shown  in  Section  2, 


1  .  »  [1  +  -  2  f  co.  ti1'2  -J0  «=«>•  ♦> 


(3-18) 


For  p£r#as  will  occur  in  the  case  of  points  lying  on  the  co-geoid  of 
the  earth, it  is  not  clear  that  the  development  which  follows  is  valid.  It  has 
usually  been  assumed  that  any  error  would  be  of  the  order  of  the  flattening  of 
the  earth,  but  recently  Molodenski  et  al  (1962)  have  shown  the  error  may  be  much 
greater. 

In  Section  2  it  is  shown  that  the  perturbing  potential  at  P  can  be 
written  as  a  series  of  spherical  harmonics  to  give 


Tp 


k  g  Yi 

k  1S2  nr+i 


(3-19) 

where  k  =  universal  gravitational  constant  and 
Yj_  =  spherical  harmonic  of  degree  i  . 

Going  to  the  spherical  approximation  in  (3-17)  i.e.,  assuming  that  we 

-V  -  d 

can  replace  by  jjrr  ,  with  sufficient  accuracy  gives 


d 

Zn 


(^)  "  3T7  [nio  fs+I  Fn  (cot 

,  ,  tl-  1 


(3-20a) 


or 


TH  (7)  “  -  n-0  ~~n-H  Pn  <co»  t) 
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( 3-20b ) 


Let  the  value  of  T  at  any  point  P  on  S  be  given  by 

•  y  i 

Tp  -  k  i^2  Tm 

r 

Again  using  ~  gives 

,Sl\_  /9T  >  _  V  f  <i+l)Yt 

W  •  (37t)  k  i *2  -~TT TT 

r  (3-20c) 

If  we  now  substitute  equations  (3-18),  (3-19) >  (3-20a),  and  (3-20c) 
into  (3-17)  ™  8**Tf  -  j-  S,  I  „10  jin  P"  <c°*  **  ill  ~77^  * 

A  pfi,  Jo  ^  '■  <«■  *>»•  (3.21) 

If  the  surface  S  were  a  sphere  (r'  =  R),  it  could  be  shown  that 


J7  Pn  Yi  dS  «  0  if  n  *  i 

*  *  (3-22) 

Thus,  again  assuming  that  S  may  be  considered  a  sphere  with  sufficient 
accuracy  that  (3-22)  is  satisfied,  yields  the  following  result.  In  (3-2l)  only 
those  terms  where  n  =  i  will  be  different  from  zero  after  carrying  out  the  indi¬ 
cated  multiplications  and  integrations.  This  result  leads  to  the  conclusions 
that: 

A.  Since  there  are  no  Y0  or  terms  the  PQ  and  P]_  terms  in  the  sum¬ 
mations  with  respect  to  n  will  contribute  nothing  to  the  final  answer.  Thus  we 
will  drop  the  PQ  and  Pj  terms  and  begin  the  summation  at  n  =  2; 

B.  The  following  identity  may  be  used 


n«2  Pn  i?2  (i  +  11  ■  J2  (n  +  l>  Pn  il2  «  ,  , 

(3-23) 

Using  the  results  (A)  and  (B)  above  and  the  approximation  r’  *  R  we 


can  write  (3-2l)  in  the  form 

.  n-  1 


T*  -  ^  [<„S2  (Ptnil  '  Pn><k  if2  ^TT>  +  <„S2  -* 


n-  1 


Pn)(k  J 


2  pr?i)lds 


(3-2U) 
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Tp 


1  rr  r  $  (2n+l)RIU  1Pn,  p  Y1  ,j„ 

TP  ■  Xi  JJ.  n«2  - T+I - 1  [k  IE2  7rn,ds  (3.25) 

Substituting  (3-20b)  and  (3-20c)  into  (3-16)  after  allowing  the  approxi¬ 
mation  r*  =  R  gives  for  every  point  on  S 

-  V  ■  ^  *  -  &  ma 

or  R  R 


(3-26) 


(3-27) 


Ag  -  (gp  -  VQ)  -  k  J2  ■■pi+2Y1 

Now  returning  to  the  integral  of  (3-25)  and  altering  its  form  slightly,  we  get 

«  -  jn  IS.  cJi  «»  *  ■>  »■  J2  <y> 

R 

Or  making  use  of  (3-23) 

*  ■  171  IS.  <A  <5?r>  "■ 

R  (3-28) 

Substituting  from  (3-26)  irto  (3-28)  gives 

'  M.  I »?2  (^TI)  Pn  (?>n+lj  A*dS  "  n-2  “ST1  (|)"+1IsA«PndS 

(3-29) 

Substituting  (3-29)  into  (3-5)  after  letting  the  point  P  lie  on  S  so 
that  gives  for  the  undulation  of  the  co-geoid  at  the  point  P. 

*  ■  S  -  ^  j2  a* 

Q  Q  (3-30) 

Equation  (3-30)  is  Stokes'  equation  relating  the  co-geoidal  height  at 
P  to  the  integral  of  the  gravity  anomalies  and  a  function  of  position  over  the 
co-geoidal  surface  in  its  spherical  approximation. 

There  are  of  course  many  ways  of  arriving  at  equation  (3-30)  (See  for 
example  Heiskanen  and  Vening-Meinesz,  1958).  The  method  chosen  here  is  meant 
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to  show  as  clearly  as  possible  what  approximations  are  made  and  where  they  enter 
the  calculations.  To  allow  for  numerical  calculations  it  is  desirable  to  have 
the  summations  of  (3-29)  or  (3-30)  in  closed  form.  In  the  following  development, 
which  follows  closely  that  given  by  Nagy  (1962),  the  closed  form  of  (3-29)  will 
be  obtained.  The  closed  form  of  ( 3-30)  would  immediately  follow  by  letting  P  =  R 
in  (3-29)  and  substituting  into  (3-5). 

To  begin  let  us  do  the  following  to  (3-29) 


a.  Factor  out 


b.  Use  the  abbreviated  notation  X  = 


c.  Use  the  identity 


from  beneath  the  integral  sign 
R 

notation  X  =  _  beneath  the  integral 

2n  +  1  ~  .  P  3 


=  2  + 


Then  the  equation  (3-29)  becomes 


Tp  ”  3712  JJ.  4g  J2  <2  +  A>  p“  xn+2< 


But  from  formula  (3-18)  using  the  approximation  r*  =  R 


(3-31) 


r  P(1  +  X2  -  2X  cos  t) 


-  Eq  XnPn 

p  n"° 


(3-32) 


\  m  J  (p0  +  P1X  +  „l2Xn  P“>  “  5  (1  +  X  cos  *  +  J2  XnPn) 


(3-33) 


Then  from  (3-32)  and  (3-33)  we  get 


nh  X"p” 


(1  +  X  -  2X  cos  f) 


172  '  1  -  x  cos  t 


(3-34) 


Using  (3-34) 

Tp  “  -A  JT  [2X2  ( - 1 - - 

4*R  (1  +  X  -  2X  cos  ♦) 

+  nil  A  p"XD+2JdS 


-  1  -  X  cos  t) 


(3-35) 
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If  both  sides  of  (3-3*0  are  divided  by  and  then  integrated  with  re¬ 


spect  to  X  we  get 

J  J,  x"' 2  pn  dx  ■  J 


X  (1  +  X  -  2X  cob  f) 


dX  -  J  — j  dX  -  J 


COS  t 


(3-36) 


Carrying  out  the  indicated  integrations  using  standard  integral  tables  yields 


J  niz  x"*2  Pn  dx  ■  nh  in  x"‘lpn 


_ dx _ _  _  (1  -  2X  cos  ♦  +  X2)1/Z 

XZ(1  -  2X  cos  ♦  +  X2)1'2  X 

2  1/2 

-  CO*  tin  (2(1.  -  —  C?6  t  ±.X  >  +1.2  cos  ♦) 


COS  t 


cos  t  In  X 


(3-37) 


Since  the  above  are  indefinite  integrals,  there  are  undefined  constants  of 
integration  which  have  not  been  indicated.  Substituting  from  equations  (3-37) 
into  (3-36)  and  using  the  symbol  to  represent  the  combined  integration  con¬ 
stant  yields 

£  1  x*1-  1  Pn  «  -  ^  -  cos  t  In  +  Y  ~  2  COB  ^  ^  “  COS  ^  In  X  +  C 

n**2  n-  1 


(3-38) 


where  to  shorten  notation  we  have  defined  v  =  (l  -  2X  cos  f  +  X2)5  -  } 
Collecting  terms  and  remembering  that  log  a  +  log  b  =  log  ab  gives 

n£2  '5TT  X"~  1  Pn  -  -  co*  ♦  In  2  (v  +  1  .  x  cot  ♦)  +  Cj 


(3-39) 


To  evaluate  ^  set  X  =  0  in  (3-39)*  Then, 


1  _n-  1 


<ni2  nfl  X  P“>x-C 


(3-to) 


[cos  t  In  2  (1  +  v  _  X  cos  f)J, 


c o 8  ^  In  4 


(3-41) 
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To  evaluate  the  function  ,  take  the  limit  as  X  approaches  zero. 

Y 

Making  use  of  the  fact  that  :  li  J 


and  noting  that 

and  that 
ve  get 


^  (x  -  2x  co®  t  +  X2)1/2  -  1/2(1  -  2X  CO*  t  +  X2) 

(-2  co*  t  + 

lim  v  *  1 
X-M) 


lim  2^!  =  lim  .  2Lz...C08  ♦ 


x«o 


x-#o 


*  cos  f 


(3-42) 


or 


Evaluating  (3-29)  at  X  =  0  yields 

0  *  cos  t  -  co8  t  In  4  + 


C .  “  cos  t  In  4  -  c o 8  ^ 


(3-43) 


Using  (3-43)  and  again  remembering  the  lavs  for  combining  log  terms,  we  get 

for  (3-39)  «  1  n  1  1  v 

n^2  nTT  X  Pn  *  "IT  “  COS  ♦  ln  2(1  +  v  .  X  cos  t)  . 

cos  if  +  cos  ^  In  4 

or 


co  y  ^  1  V 

Z  — — — w  Pn  *  -=- - cos  t  (1  +  In 

n*2  n_  l  X 


(1  +  V  -  x  cos  *)v 
- 2 - >  (3-44) 


Using  (3-44)  and  noting  that 
equation  (3-35)  becomes 


g  X 


n+2 


,n-  1 


hS2  —577  Pn  -  X3  J,  —-j  P„ 


TP  "  ^ 8  12x2  ^  '  1  '  x  cos  ♦)  +  3X3  <2^!  .  cos  ♦(!  + 


or  collecting  terms 

IP  -  ~2  !JS  X2  [^  +  1  -  3 

4nR  8  v 


l-f1  ±  v  -  X  cos  ♦ 


(3-45) 


>>>]dS 


v  -  5X  cos  f  - 


3X  cos  *  In  i-  +  v  “  x  C08  ♦) 

T~ - — 


JAgdS 


(3-46) 


Again  note  that 


(1  +  X  -  2X  cos  *)  1/2  «  I 

P 


-  1/2 
2X) 


v 


(3-47) 
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TP 


By  substituting  this  into  (3-46)  and  remembering  that  X  =  5.  we  get 

57!  s  R[l  +  ?  -  T  -  c0*  ♦ 


P  P' 

-  H  co.  *  m  l£-±J^  *  c°»  ♦)]  Agds 

P  f  3-48) 

If  we  now  assume  that  point  P  is  also  on  the  surface  S  which  can  be 

treated  with  sufficient  accuracy  as  a  sphere  then  p  =  R.  Making  the  spherical 
assumption  we  can  relate  the  distance  r  which  is  now  simply  the  chord  between 
two  points  on  the  sphere  to  R  and  t  by  the  well  known  formula 

r  *=  2R  sin  £  ♦  (3-49) 

Using  the  two  relations  p  *  R  and  r  =  2R  sin  £  f  equation  (3-46)  becomes 

Tp  ■  57!  IT.  t-.£n-172-|  +  1  -  6  *ln  1/2*  -  5  cos  t  -  3  cos  t  In 

(3-50) 


(tin  1/2*  +  9in  l/2,>)1  d8ds 


Or  in  shortened  form 


Tp  ■  571  S(*>  AgdS  (3-51) 

where  the  symbol  S(  t  )  represents  the  quantity  in  brackets  in  equation  (3-50). 

Then  substituting  (3-5l)  into  (3-5)  gives 

Np  ■  ??  "  57!y-  JJ*.  S(*>  AgdS 

Q  w  (3-52) 

This  is  the  normal  Stokes'  equation  for  co-geoidal  undulations.  Writing  equa¬ 
tion  (3-4Q)  in  shortened  forn}  gives  for  a  point  not  on  S 

*P  “  57!v“  JI  s  <p»*>  igds 

5«R7q  (3.53) 

where  S  (p,  t )  is  used  to  represent  the  quantity  in  brackets  in  equation  (3-48). 

In  classical  theory  one  carries  out  numerical  integration  of  (3-52) 
over  the  surface  of  the  earth  treated  as  a  sphere.  The  solutions  are  set  up  to 
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give  the  geometrical  quantities  represented  by  S  (f)  in  predetermined  form  far 
multiplication  by  the  Ag  values. 


Deflections  of  the  Vertical  in  Classical  Theory 

The  deflection  of  the  vertical  in  a  particular  direction  at  a  point  P 
on  the  bounding  equipotential  surface  is  the  angle  between  the  normal  to  the 
bounding  equipotential  at  the  point  and  the  normal  to  the  potential  field  U  at 
the  point.  This  will,  of  course,  be  the  angle  between  the  two  equipotential  sur¬ 
faces  at  the  point.  It  is  usual  to  express  the  deflection  of  the  vertical  in 
two  components,  an  east-west  component  and  a  north-south  cocponent.  It  i3  con¬ 
ventional  to  call  the  north-south  component  g,  and  the  east -west  component  \ 

If  N  is  computed  from  Stokes'  formula,  then  the  deflection  which  is  the 
rate  of  change  of  N  in  a  direction  say  1  (angular  direction  Ip)  is 


deflection  .  ^  N  -  i  (Ij)  ~  ~  <TP> 


(3-5*0 


Some  confusion  can  arise  as  to  the  sign  of  the  deflections.  The  sign 
is,  of  course,  a  convention  since  the  deflection  is  simply  a  difference.  The  normal 
convention  is  to  define  the  meridional  deflection  g  as  the  astronomic  latitude 


minus  the  gravimetric  latitude  and  the  parallel  deflection  T\  as  the  astronomic 

longitude  minus  the  gravimetric  longitude. 

i.e.  g  3  •  *st  -  9  grav. 

T]  =  X  ast  -  X  grav. 

Let  us  consider  the  meridional  deflection  g  and  the  picture  below 

qtoid. 
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but  dS  =  R  sin  fdf  dA 

5 


cot  or  tin  f  dAdf 


“  ^  ^"[jf  3(^]4* 

11  "  K  ^“[yf  8<t)]4g  ,ln  *  ,ln  ♦  <«*<»♦ 


(3-59) 


(3-6o) 


where  dA  is  an  incremental  angle  corresponding  to  a  longitude  angle  using  the 
computation  point  as  a  pole. 

Let  us  evaluate  h  sit) 


yj  s(t) - z 


_.<r°*  i/ -1  .  3  cos  1/2*  +  5  sin  *  +  3  tin  *  In  («in  1/2*  + 
sin  1/2* 

2  .mix.  3  COB  ♦  cos  l/2t  (1/2  +  tin  l/2f) 
tin  l/2t)  7ln  1/2^  (1  +  tin  1/ZfJ 

If  we  now  multiply  through  by  sin  t  and  use  the  identities  sin  t  =  2  sin  \  ♦ 

cos  £  t  and  cos  ♦  «  1  -  2  sin^  t/2  we  get 

.in  *  yr  S(t)  -  .  6  .In  |  (1  -  Bin2  + 

5*  Bin  1/2*  2  2 

20  Bin2  |  (1  -  Bin2  ^)+  3  ,in2  f  In  (»in  1/ 2t  +  Bin2  1/2*)  - 

6(1  -  2Bin2  |)(1  -  Bin2  |)(l/2  +  Bln  1/2*) 


Bin  *  8  (*)  -  -  esc  |  -  3+3  sin2  *  In  (sin  1/2*  +  sin2  1/2*) 

-  8  sin  |+32  Bin2  ^  +  12  *ln3  \  -  32  Bin4  |  (3-6l) 

Then  if  we  let  F'(*)  =  £  sin^S(*)  we  get 


s 


r  r*  '■(♦>  '«• « <»Ad* 


7*Yq  0  0 

p2n 


11  -  f  J2*  F'(*)  Ag  sin  a  dAd* 


Q  0  0 


(3-&) 


We  might  note  here  that  for  small  values  off cosf«l,  sin  and  thus 
we  can  use  the  approximation  F  m  “  1/2  (esc  ^  +  3) 


,ln  ^  Tf  ®W)  "  "  c,c  -  3 


(3-63) 


Numerical  methods,  of  carrying  out  the  integrations  of  equations  (3-62) 
have  been  developed  by  a  number  of  workers.  These  numerical  methods  can  be,  to 
some  extent,  carried  over  to  the  method  of  solution  adopted  for  use  in  the  Test 
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and  Application  Phenes  of  this  contract.  Ihe  overall  numerical  method  finally 
adopted  is  discussed  in  Section  7* 
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SECTION  k 

The  New  Geodetic  Theory 

From  the  classical  viewpoint  described  in  the  previous  section,  the 
necessary  steps  in  surface  level  deflection  computations  are  (l)  correct  transfer 
of  material  from  above  to  below  sea  level,  (2)  "reduction"  of  surface  gravity 
anomalies  to  geoid  and  then  co-geoid  level  anomalies,  (3)  computation  of  deflec¬ 
tions  at  co-geoid  level,  and  (4)  computation  of  surface  level  deflections  from 
the  co-geoid  level  deflections. 

The  above  procedure  appears  to  result  in  considerable  inaccuracy 
because  of  lack  of  knowledge  concerning  internal  densities  and  vertical  gravity 
gradients.  Because  of  these  inaccuracies  and  the  difficulty  involved  in  assessing 
the  errors  introduced  from  approximations  made  in  the  mathematical  development  of 
the  theory,  considerable  argument  exists  as  to  the  "correct"  reduction  to  be  used 
and  the  type  of  anomaly  which  is  theoretically  more  accurate.  It  has  not  been 
possible  until  recently,  because  of  the  way  in  which  the  theory  was  developed, 
to  judge  accurately  which  reduction  method  is  best. 

The  so-called  new  theory  of  gravimetric  geodesy  developed  below  is 
therefore  important  for  two  reasons. 

1. )  It  makes  possible  an  evaluation  of  the  accuracy  of  the  results 

obtained  using  different  types  of  gravity  anomalies  with  the 
classical  theory. 

2. )  It  provides  a  mathematical  theory  capable  of  extension  to  any 

order  of  accuracy  desired. 

The  following  development  of  the  new  theory  follows  very  closely  that 
of  Molodenski  et  al  (i960).  First  we  will  develop  an  equation  analogous  to  Brun’s 
equation  for  the  co-geoidal  height  in  classical  theory. 


» 
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Let  W  (0,  X,  r)  =  potential  field  of  the  actual  earth  Including  the 

the  rotation  potential 

U  (9,  X,  r)  *  reference  potential  field  having  the  same  rotational 
potential  as  the  earth's  field 

WQ  =  value  of  sea  level  equipotential  of  the  actual  earth. 

Now  assume  we  pair  the  equipotential  surfaces  of  the  actual  earth  and 

the  reference  field  in  some  way  so  that 

UQ  =  the  equipotential  surface  of  the  reference  potential 

field  which  is  paired  with  equipotential  W0  of  the 
actual  earth. 

The  most  advantageous  way  to  pair  the  surfaces ,  is  to  pair  the  equipotentials 
which  have  the  same  value  of  potential.  If  this  is  done  we  have 

u0  =  w0  (4-1) 

We  cannot  do  this  though,  without  knowing  WQ.  As  is  the  case  with  the  classical 
theory,  the  assumption  is  made  that  WQ  is  known  although  it  is  not  in  fact  known 
exactly. 

Consider  some  point  P  (0*,  X*,  H)  on  the  earth's  surface  where  H  is  the 
height  of  the  point  P  above  the  surface  U  =  UQ  whose  shape  is  known.  Let  us  de¬ 
fine  the  perturbing  potential,  T,  at  any  point  to  be  the  difference  between  the 
actual  earth's  potential  at  the  point  and  the  reference  potential  at  the  point. 
Thus 

T  <9*,x*,h)  -  W  (S*,X*,H)  -  u  <e*,H)  (4-2) 

p  P  r 

Note  that  this  equation  assumes  that  U  has  been  chosen  so  as  not  to  be  a  function 
of  the  coordinate  X  .  This  is  the  case  for  the  international  formula  potential. 
Then 

w  (9*,X*,H)  -  W.  -  Jgdh 

p  (4-3) 
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where  g  is  the  actual  value  of  gravity  along  the  path  of  integration,  and 
integration  takes  place  along  the  earth's  surface. 

Let  us  now  define  the  equipotential  surface  Uq  (0fh)  which  is  paired 
to  Wp  by  the  requirement  that 


Wo  -  U0  =  Wp  -  Uq  (k-h) 

Then  we  have 

UQ  “  u*  -  /Ydh  (4-5) 


where  J’ydh  is  defined  by 

J*  gdh  *  J  Ydh  -  hy 


(4-6) 


y  =  mean  value  of  y  along  line  of  integration. 

In  principal  h  can  be  rigorously  determined  by  gravity  and  leveling  on 
the  earth's  surface.  In  practice  h  is  not  so  rigorously  determined  due  to  the 
lack  of  gravity  information  along  the  survey  lines.  However,  the  difference  be¬ 
tween  normal  map  elevations  and  the  precise  quantities  is  so  negligible  as  to  be 
unlikely  to  cause  any  significant  error.  From  this  point  on  the  symbol  h  will 
be  understood  to  represent  normal  orthometric  heights  as  read  from  maps. 

Now  by  substituting  for  Wp  from  (4-4)  we  can  write  (4-2)  in  the  form 

Tp(9*.  X*,H)  -  W.  -  U.  +  U  (6,h)  .  U  (8*,H) 

(4-7) 

As  the  differences 

9*  -  0  ■  &9 
X*  -  X  -  AX 


are  small  quantities,  their  squares  and  derivatives  can  be  ignored.  If  we 
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expand  U  about  the  point  Q  in  a  Taylor  series  keeping  only  the  first  terras, 
obtain 


up(e*,H)  -  uq(9,h)  +  (|£)q  h  +  (yg)p  A6 


,au, 


we 


(4-9) 


Now  the  error  in  A0  is  normally  less  than  1*  so  that  the  last  term  in  (4-9)  may 
be  ignored.  Substituting  (4-9)  into  (4-7)  and  ignoring  the  last  term  in  (4-9) 

we  get  *  „ 

Tp  -  W.  .  U.  +  UQ  -  0Q  -  (|")QH 


or 

Tp  -  W.  -  U.  .  (|«)QK 

(4-10) 

But  to  a  very  close  approximation  =  gravity  due  to 

normal  potential  field  at  Q.  This  gives,  after  rearranging  terms  and  substituting 
into  (4-10), 


.  Tp  W„  -  U, 


(4-11) 


If  we  have  WQ  ■  Uo>  we  obtain 


(4-12) 


which  is  completely  analgous  to  Brun's  equation.  Here,  however,  Tp  is  the  de¬ 
viation  potential  at  the  earth's  surface,  is  normal  gravity  at  a  point  distance 
N  above  or  below  the  earth’s  surface.  N  is  thus  a  correction  term  to  be  added 
(with  correct  sign)  to  h  to  give  true  location  above  the  equipotential  U  *=  Uq  of 
the  surface  point  P.  N  is  usually  called  the  height  anomaly  in  this  case. 

The  assumption  that  WQ  =UQ  made  in  deriving  (4-12)  from  (4-11)  has  been 
examined  in  detail  by  Mo lode ns ki  et  al  (i960). 
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The  solution  of  equation  (4-12)  requires  knowledge  of  T  =  W  -  U  .  To 

P  P  P 

obtain  Tp  let  us  first  establish  the  boundary  condition  which  T  must  satisfy  at 
the  earth's  surface.  Since  W  =  U  +  T  we  can  write 


(SVp 


Wp  + 


v=  direction  normal  to  reference  equipotential 

Now  with  an  error  of  the  order  of  the  deflection  of  the  vertical,  we  have 
(^J)p*gp=  observed  gravity  at  point  P.  Expanding- about  Q  in  a  Taylor  series 
and  keeping  only  the  first  term  we  get 


Substituting  this  gives 


8p  ■  VQ  '  M4^)Q  +  (^}p 

or, using  Brun*s  formula,  we  get 

(8p  ■  v*  +  4^? 


,dT, 


(U-13) 


Thus  the  determination  of  T  has  been  reduced  to  the  solution  of  the  third 
boundary  value  problem  of  potential  theory  at  the  earth* s  surface.  Ms  problem 
can  be  stated: 

Find  the  function  T  which: 

1.  )  At  every  point  on  the  physical  surface  of  the  earth  S  satisfies 

the  condition  (4-13) 

2. )  Is  a  harmonic  function  outside  S 

3. )  At  infinity  is  regular  or  satisfies  the  stronger  condition 

piiS.-  (T,2)  =  o 

To  find  the  function  T  it  is  possible  to  proceed  in  a  number  of  ways. 

Some  authors,  e.g.  Hirvonen  (i960),  have  proceeded  much  along  classical  lines  by 
using  spherical  harmonics  and  arriving  at  the  generalised  Stokes  and  Vening  Meinesz 
theorems.  Molodenski  et  al  (1962)  have  examined  these  methods  and  point  out  that 
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the  error  is  greater  than  the  order  of  the  flattening  of  the  earth  as  has  often 
been  supposed.  They  also  stress  that  the  error  in  deflection  computation  is 
much  greater  than  in  geoid  undulation  computation. 

Another  method  of  procedure  and  the  one  favored  by  Molodenski  involves 
the  introduction  of  a  surface  density  coating.  The  use  of  density  coatings  dates 
back  to  Helmert  and  is  based  on  that  part  of  potential  theory  which  permits  the 
potential  of  a  volume  mass  to  be  represented  external  to  it  by  the  potential  of  a 
surface  mass  distribution  on  the  surface  of  the  volume  mass. 

Thus  to  solve  equation  (4-13)  we  introduce  the  auxiliary  density  function <p 
representing  this  surface  density  and  defined  by 


T 


P 


(4-14) 


where  rp  is  the  distance  from  the  point  P  where  T  is  being  computed  to  the  variable 
point  on  the  surface  S. 

To  substitute  into  (4-13)  we  also  need  (|*)p«  As  the  point  P  lies  on 
the  surface  S,  the  derivative  contains  an  extra  term  due  to  the  density  layer 
causing  a  discontinuity  of  the  derivative  at  the  surface.  Thus  we  arrive  at 


+Op  ■  +  2n»p c°*  *  +  d"iP 


(4-15) 


where:  +2tto  cos  Of  is  the  term  appearing  due  to  the  discontinuity  of  the 

P derivative  at  the  density  layer  S 

and  Of  is  the  angle  between  the  normal  to  the  direction  in  which  the 

derivative  is  taken  (-Lto  the  reference  field)  and  the  normal 
to  the  surface  S. 

If  we  can  derive  a  tp  and  satisfy  equation  (4-13)  by  the  T  determined 
from  it  using  (4-14)  then  this  T  will  obviously  satisfy  boundary  conditions  2.) 
and  3»)  i.e.,  T  will  be  harmonic  outside  S  and  regular  at  infinity.  This  results 
from  the  fact  that  a  T  so  obtained  will  be  the  potential  of  an  actual  density 


distribution  and  must  therefore  behave  as  desired  outside  S. 
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To  do  this  we  can  substitute  (4-14)  and  (4-15)  into  equation  (4-13) 


which  allows  us  to  write  this  equation  in  terms  of  <p  .  Substituting  we  get 

2tt®  cos  a.  +  f  fP-  dal  -  i  f  L  ds  -  (g  -  y  ) 

+  l5vJs  rp  Jp  yqWqJ«  rp  VKP  V  (1|  l6) 

where  r^  is  the  distance  from  the  point  P  to  the  variable  point  on  S. 

Since  the  surface  S  is  unknown  the  equation  (4-i6)  is  unsolvable  in 
the  form  given.  However,  a  sufficiently  reliable  approximate  solution  can  be 
obtained  by  carrying  out  the  operations  indicated  on  the  l.h.s.  of  (4-l6)  on 
the  surface  S  which  is  the  first  approximation  to  the  earth  obtained  by  adding 
the  measured  heights  h  to  the  international  ellipsoid.  Then  we  get 

'2%  — » +  H£-q  d«  -  \  di-  <*p  -  V 


and  the  approximation  of  Tp  is  given  by 


Tp  -  J-  I"  di 

8  rQ 


(4-18) 


If  we  designate  the  radius  vector  of  the  point  Q  by  ^  ,  the  radius  vector  of 
the  element  ds  by  p,  and  the  angle  between  the  two  radius  vectors  by  f,we  have 


Pq  +  P2  -  2ppq  c0*  ♦ 


If  we  assume  that  the  reference  potential  is  spherical  we  have 

A  dYx  -  z1 


(4-19) 


(Y  Sv}Q  ""  (Y  ^p)Q  "  + 


and 


2 

PQ 

-  3  ,1 


&  (VJQ  (^)lQ 


(4-20) 

-3/2 


7  2  1/2  22 

or  [(p*  +  pQ  -  2pp  Q  cos  ♦)"  1  «  -  1/2  (P  +PQ  -  2PPq  cos  f) 

Q  2  2  2 

2p  -  2p  p  cos  t  P  -  P  o  I 

( 2 p  -  2p  cos  t)  -  -  1/2  [ - 2 - 3 - ]  -  - -  - 

Q  P-r  -  2prtr  „  2p  r 


Q  Q 


Q  Q 


Q  Q 

(4-21) 
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Then  substituting  (4-20)  and  (4-21)  into  (4-17)  we  get 


2*  co.  or  -  (gp  -  Yq)  +  ^  J  i  d8  +  ^  J* - -T-1 - 


(4-22) 

What  ve  have  assumed  here  has  nothing  to  do  with  the  surface  of 
integration.  We  have  merely  assumed  that  the  equipotentials  of  U  are  sufficiently 
near  spherical  that  equations  (4-20)  and  (4-21)  have  sufficient  accuracy  for 
practical  purposes. 

Equation  (4-22)  is  a  linear  integral  equation  in  9.  The  method  of 
obtaining  a  solution  to  the  problem  is  as  follows . 

1. )  Solve  (4-22)  for  9. 

2. )  Use  9  to  obtain  T  from  (4-lQ) 

3. )  Use  T  to  obtain  height  differences  (N)  from  (4-12) 

4. )  To  obtain  deflections  of  the  vertical  ve  would  need  the  formula 

S3  “  7Z  J  ?  ds  -  2n*  coi  <">“> 

where  m  is  a  direction  perpendicular  to  the  radius  vector  of  a  point. 

As  ve  have  made  the  spherical  assumption  for  the  reference  potential,  let 
us  denote  the  spherical  reference  surface  as  having  a  radius  R  then  we  have 


R  +  h 


P  -  R  +  h 


(4-23) 


with  an  error  of  the  order  of  the  flattening  of  the  earth. 

Then  with  only  an  error  of  the  order  ^  which  is  permissible  since  we 
are  retaining  only  accuracies  of  the  order  of  the  flattening  we  have 


p2  .  P2q  -  2hR  -  2hqR  +  h2  -  h2q 


(4-24) 


or 


.  p  Q  «  2R  (h  -  hQ) 


and 


(4-25) 
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This  gives 

2*»q  cos  or  -  (gp  -  Yp>  +  «  J.  f  ds  +  J.  T5®  *ds 

Q  r  Q  (4-26) 

From  this  point  one  could  proceed  directly  to  the  iterative  solution  of 
Molodenski,  but  it  is  perhaps  instructive  to  see  how  the  first  approximation 
which  is  the  one  most  often  used  arises.  The  following  development  is  modeled 
after  Arnold  (i960). 

If  we  examine  the  method  of  solution  of  the  differential  equation  (4-13) 
to  arrive  at  Stokes'  equation, we  see  that  what  is  assumed  is  that  the  integra¬ 
tion  is  carried  out  over  a  sphere  and  that  both  the  computation  point  and  the 
variable  points  lie  on  this  sphere  (  h  =  h^  =  0).  Thus  the  Stokes  equation 
approximation  amounts  to  solving  (4-26)  in  the  simplified  form 


2*  <Pp  -  (8-V)  +  A  Ac  ^  dX 


(4-27) 


where  X  denotes  integration  over  the  spherical  surface  and 

dX  ■  sin  fdfdA  . 

We  have  already  shown  that  starting  with  the  differential  equation  form 
of  (4-27)  we  arrive  at  the  solution 


1  *  ?Tr  lx  (8P  ■  V  s  <cos  ♦>  dx 

Thus  we  get  from  (4-14) 

T  “  J*.  ds  -  -T*  (gp  -  V  s  <c°s  ♦)  dx 

where  S  (cos  f)  is  the  Stokes  function. 

Then  from  (4-27)  we  have  as  a  first  approximation  to  9 


(4-28) 


(4-29) 


<p, 


(8p  -  V  1  3 

TOT 


Jlc  *8p  ■  v  s  *cos  ^  dx 


(4-30) 
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Now  to  arrive  at  a  secondary  approximation  we  go  back  to  equation  (4-26). 
The  spherical  assumption  is  retained  for  integration  of  the  first  integral  on 
the  r.h.s.  of  (4-26).  In  the  second  integral,  the  approximation  h  =  hQ  =  0  is  not 
made  but  it  is  assumed  to  be  accurate  enough  to  replace  cp  in  the  integral  by  <p^ 
from  (4-30).  This  gives  the  equation 


2*<p  -  (8  -  V)  +  7l-  L  «Pt  dX  +  Ti  J“x  f 


dX 


(4- 31a) 


ill— h  q 

If  we  assume  that  within  the  integral  J- — 3— ■  cp^  dXit  is  accurate  enough 
to  replace  Vy  (S2^Q)  rather  than  the  more  accurate  expression  from  (4-30), 


we  get 


2**  -  (gp  -yq>  +  ^  fx  (gp  .  Yq)  dx  +  ±  Sx  ^dx 


(4-31b) 


By  analogy  with  the  solution  of  (4-27)  we  see  that  if  we  consider  our 

1  ii  h-  h 

'’anomaly"  to  be  (gp  -  Yq)  +  - g  (gp  -  Yg)dX  in  (4-31) 

as  compared  to  simply  (g-y)  in  (4-27)  we  should^arrive  at  the  same  type  of 
solution  as  (4-28)  i.e., 

T  -  ^  fx  I(8p  -  YQ)  +  AS*1S  (cos  ♦)  dX 

(4-32) 


where 


Ag* 


-1-  f 

2 it  Jx 


Yq)  dX 


(4-33) 


Equation  (4-32)  is  the  equation  for  T  from  the  new  theory  most  commonly 
used  in  arriving  at  corrected  formulae  for  deflections  of  the  vertical.  The 
method  of  arriving  at  equation  (4-32)  as  a  second  approximation  may  not  seem  as 
rigorous  as  desirable.  However,  one  can  arrive  at  this  same  result  in  a  more 
rigorous  manner  by  proceeding  as  Molodenski  et  al  (i960)  have  done.  The  para¬ 
graphs  below  outline  their  procedure. 


We  shall  start  from  equation  (4-22)  which  we  repeat  here  as 


2*<p  CO.  •  -  <s  -  vQ)  +  J  A  ds  +  TT  “ 


ft2  A  2 
P  "po 


(4-34) 


where  we  have  removed  the  bars  to  denote  another  surface.  Therefore,  we  must 
remember  that  in  what  follows  the  surface  S  in  (4-34)  is  the  first  approximation 
to  the  actual  earth's  surface  which  is  obtained  by  adding  the  normal  height  to 
the  international  ellipsoid. 

Let  S  of  equation  (4-34)  be  the  first  approximation  to  the  physical  surface 
of  the  earth.  Suppose  we  have  another  surface  S  which  is  related  to  S  by  the 
following  transformation.  The  angular  coordinates  of  the  corresponding  point  on 

S  are  the  same  as  the  angular  coordinates  of  a  point  on  S.  The  radius  vectors 

of  corresponding  points  are  connected  by  the  equation 

p  -  R  +  k  (p  -  R)  -  R  +  kH 

(4-35) 

where:  P  =  radius  vector  to  transformed  surface  S 

k  *  some  constant  coefficient 
p  =  radius  vector  to  surface  8 
R  =  radius  of  mean  sphere  of  earth 
H  =  normal  height 

If  we  let  k  =  l  this  would  correspond  to  transforming  3  into  itself.  If 
k=0  the  surface  S  is  a  spherical  surface  of  radius  R. 

Now  let  us  go  back  to  equation  (4-34)  which  gives  us  an  equation  at  the 


earth's  surface 


2icf0  co*  a  "  (8p  '  V  +  TJ~  J.  r  d8’  +  TT"  J*. 


Let  us  define  a  new  surface  density  function  by  the  relation 


fl2  2 
P  -po 


<p  dS' 

(4-36) 


f  J.C  «  or  f 


X  ±2  eo.  a 


X 


(4-37) 
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R2  X 


Substituting  into  (4-34)  we  get  2jt  X  cos  ar«(g-Y)  +  jp-  J8  -7  * 


2  2 

P  -P. 


®  dS  +  J*g  - 3 - X  — ^  cosefdS 

of  fl 


If  we  now  revert  to  the  incremental  solid  angle  by  the  relation 
dw  =  — Sxj~  ^  assuming  S  is  near  enough  a  sphere  to  eliminate  integration  with 


respect  to  p  and  divide  both  sides  by  R  we  get 


2*  -*7  cos2  .  -  -  ‘  -  (g  -  Y)  +  jI~  /  S  - w  +  ,i-  J 

P  R  o  or 

°  2 
If  we  now  multiply  both  sides  by  pQ  we  get 

2  3p  2  2 

2k  X  cos2  Of  -  (8  -  V)  +  -y-  f  *  dv  +  f  - 3—  xdw 


(4-38) 


In  this  expression  remember  that  dw  is  the  element  of  solid  angle  subtended  at 
the  origin  by  a  surface  element  dS. 

Now  let  us  consider  a  surface  S  related  to  the  surface  Sfby  equation 
(4-35)  with  l^k^O.  If  the  anomalies  (g  -y)  had  been  observed  on  this  surface 
we  could  write  an  equation  analogous  to  (4-38)  in  terras  of  "barred"  quantities 


related  to  this  surface.  Thus  we  would  have 


2j(X  cos2  a 


P-y  (8  -  Y)  +  \  PQ  J  *  dw  +  \  pQ  / 

R  r 


n2  n  2 

P  -P 

_-=3£_  xdw 


(4-39) 


Equation  (4-38)  will  have  a  solution  provided  (g  -  y)  satisfies  the 
condition  J'(g-Y)Yidw=0.  When  going  from  S  to  S  the  element  dw  does  not 
change  and  the  anomalies  (g  -  Y)  are  the  same.  So  if  (4-38)  has  a  solution  so 
does  (4-39)*  The  perturbing  potential  at  the  surface  S  for  the  situation 


given  by  (4-39)  becomes 


R  J  l 


(k- ho) 


To  recapitulate  what  we  have  shown  is  that  if  the  gravity  anomalies  (g  -  y) 
belonged  to  the  surface  S  related  to  s'by  (4-35)  rather  than  to  S,  the  external 
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perturbing  potential  would  be  given  by  (4-40)  where  X  was  a  solution  of  (4-39) 
We  can  show  by  geometric  reasoning  that 


f2  -  rQ2  (1  +  k  +  k2  ™2)  +  k2  (H-H0)2 


(4-4l) 


where  r,  is  the  distance  between  the  point  under  consideration  on  S  and  the 


moving  point  on  S,  and 

rQ  *  2R  sin  ^ 


(f  =  angle  between  p  and  p  ) 

0 


We  note  that  by  geometric  reasoning  it  can  also  be  shown  that 


tan2  *  ■  (H>2  +  <tf>2 " 1(2  tani  " 


(4-42) 


Now  if  we  expand  X  in  powers  of  k  i.e.,  replace  X  by  a  series,  then 


|  X  k' 
o  n 


(4-43) 


where  the  functions  Xn  are  independent  of  k. 


Since  the  dependence  on  k  of  the  functions  p,  r,  cos  a  are  given  by  (4-35) 
and  (4-42),  we  may  now  substitute  into  (4-39)  and  expand  both  sides  of  the 
resultant  equation  in  powers  of  k.  The  resulting  series  appearing  on  the  two  sides 
of  the  equation  must  then  be  identical  regardless  of  the  value  of  k  (l>kX>);  thus 
the  coefficients  of  each  kn  term  in  the  left  and  right  hand  sides  of  the  equation 
must  be  equal.  We  thus  obtain  an  infinite  system  of  integral  equations  for  which 


one  obtains  the  functions  XQ,  X-^,  X^,  etc.  successively  for  as  many  terms  as  are 
desired.  The  results  then  hold  for  any  surface  S  where  l>k>0  .  In  particular 


the  results  will  hold  for  bwl,  i.e.,  for  the  first  approximation  to  the  physical 
surface  of  the  earth. 

Now  we  will  solve  (4-39)  after  making  the  expansions.  First  note  that 
(4-39)  has  a  relative  error  of  the  order  of  the  flattening  of  the  earth  since  we 
have  taken  the  reference  figure  to  be  a  sphere  rather  than  an  ellipse  of  revolution 
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in  developing  the  formula.  Thus  ve  will  simplify  the  calculations  by  allowing 
errors  of  equal  magnitude,  i.e.,  errors  of  the  order  Jj  .  This  will  allow  us  to  make 


the  approximations 


?o  "  R 


Then  we  get 


r2  -  rf  +  k2  (H- H0 ) 2 


(4-44) 


cos  5  -  (X  +  k2  tan2  i 


(4-45) 


2«*  co*2  5  -  (8  -  V)  +  \  R  J“-§-  dv  +  R2;  Xdw  {k.k6) 

r  r 

Now  expanding  in  series  ve  get 

2n  (X  +  kx,  +  k2x,  + - )(1  -  k2  tan2  or  +  k4  tan4  or - )  -  (g  -  y) 


*  J*  ~  (1  - 

L  r\ 


1  .2.  2  ,  3  4.4 

-*  k  h  +  k  h 


)  (Xn  +  kX  1  +  k  + - ) dw 


+  R  J  —  3  <H-H0)(1  -  \  kV  +  ^  kV  +  ----  )(X0  +  kXx  +  k2X2  + - )dw 

j.© 

where  h  =  (H-Ho)/r0 

Equating  the  coefficients  of  each  power  of  k  we  get 

2*Xl  -  |  R  J  |i  dw  +  R2  J  dw 


2rtXo  -  (g  -  V)  +  \  R  J*  f8  dw 


(4-48) 


etc.  with  similar  equations  for  higher  order  terms.  We  note  that  all  of  these 
equations  have  a  solution  exactly  analogous  to  the  solution  of  (4-27)  with  the 


"anomaly"  having  different  forms.  Thus  in  an  analogous  manner  we  are  lead  to 
solutions  l  f“  dw  ‘  *7R  J*  (g  -  Y)  S(cos  ♦>  dw  „  .  . 

.  (4-49) 


-  +  — --7  I  («  -  y)[s(co«  ♦)]. 


(4-49) 


Then  using  this  ve  get  ;  ^  ^  J  ,  R2  J  2^  f# 

[V  J  S,2°  Xo  dJ 

X  -  _ _ _ _ *■  +  — ““7  J  l  R  J*  ““■§  x0dw J  [  S  (co 8  ♦  )  ]dw 


J  lR  /  ““ — §  x0dwl  [  S  (cos  t)]i 


(4-50) 


Molodenski  et  al  (i960)  write  the  general  equation  as 


whose  solution  is 


with 


-  „  .  -Hr  I  --  d«  ”  °n 

2«Xn  -  J  rQ 

/  dw  -  J*  Gn  S(coa  f)  dw 
Xn  -  3  +  ---2  J  Gn  S(coe  ♦  )  dw 

(4*) 


G]  =  r2  r  x  dw 

1  «  J  o 


(4-51) 


(4-52) 


(4-53) 


Go  ”  (8  -  Y) 

Then  to  get  values  of  T  we  return  to 

T  =  R2  J*  |  dw 

and  carry  out  the  series  expansions  we  arrive  at 

-  n2  r  1  /1  1  I2*2  ^  3  ,4  ,4  ,  x 

T  -  R  J  ~  (1  -  -2  4  h  +  3  4  h  + - ) 

<Xo  +  kXl  +  42X2  + - )dw  (4-54) 

where  again  remember  that  h  =  **  -  *4) 

r0 

On  collecting  powers  of  k 

x  -  R2  /  |u  dw  +  kR2  J  ~  dw  +  k2R2  (J  -  J  J  -(5-^l->2  Xo  dw] 

r°  (4-55) 


+ - 

where 

T 


£  k"T 
o  n 


etc . 


P  G  S(cos  ♦)  dw 
j  o 

T1  *  J*  G1  s(C0S  ♦)  dw 


R 
o 

R 


(g  -  Y) 
R 2 


X  dw 

o 


To  get  the  value  on  the  first  approximation  surface  we  would  of  course 
set  k*l.  If  we  set  k  =  0  we  have  the  result  which  would  have  been  obtained  by 
letting  all  the  anomalies  be  referenced  to  a  sphere. 

We  have  in  the  development  assumed  that  Uq  =  Wq  and  thus  do  not  have  the 
constant  term  l/2  in  the  integrals  as  Molodenski  et  al  (i960)  do. 
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Molodenski  et  al  (1962)  point  out  that  the  equations  for  the  potential 

give  only  the  potential  at  the  earth's  surface.  Thus  we  know  only  the  quantity 

Ts  =  potential  at  the  earth's  surface.  In  determining  the  deflections  we  wish  to 

compute  the  rate  of  change  of  T  in  a  circular  direction  1.  The  deflection  In  this 

direction  can  be  expressed  as 

0  -1  dT  -1  dT 

Yq  ‘  y?  ** 

(4-56) 


where  dl  =  Rd<p 


dT 

To  get  at  a  point  in  terms  of  Tg 


we  note  that 


d  Ta  ST  ,  9T  9h 
59  59  51  59 


Thus 


dT  aT.  st  3h 
5 9  "  &9  3h  ^9 


dT 

To  evaluate  we  use  equation  (4-13)  which  we  approximate  by 
dH 

up  V  ~  -  Sr 

dT  dT 

Using  the  further  approximation  m  gg*  this  gives 

51  Up  V  +  ~T 


Substituting  into  (4-58)  this  gives 


9T  ST,  f.  v  i  +  2Tpl  Sh 

59  "  59  ,(8p '  V  +  nr1  39 


This  gives 


q  “ 1  1  ST.  .  f  .  .  2Tp,  SH  1 

6  "  y~  jl5?*  ((8p  -  V  +-TT1  39  1 


P  '<T 


59  R 


(4-57) 


(4-58) 


(4-59) 


(4-60) 


(4-61) 


(4-62) 
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Or  substituting  for  Ts  from  (4-32) 

•  ■  3-  Irarf,  «».  -  v„)  .  ».*)  H  ||  «  t  t<«„  -  v.)  ♦  3*1  i  IS 


Yq 


* 

(4-63) 


This  is  the  basic  deflection  formula  as  developed  in  the  new  theory.  It 
consists  of  the  Vening  Meinesz  formula  plus  two  correction  terms.  The  first  cor¬ 
rection  consists  of  the  Ag*  term, which  is  a  correction  for  the  fact  that  during 

integration  a  sphere  is  used  instead  of  the  correct  surface  of  the  earth  in 
dH 

determining  T.  The  term  takes  into  account  the  fact  that  in  deflection  compu¬ 
tation  we  need  to  get  derivatives  in  directions  other  than  along  the  earth's 


surface . 
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SECTION  5 

Theoretical  Aspects  of  Anomaly  Selection 

In  this  section  we  shall  examine  the  theoretical  aspects  of  anomaly 
selection.  That  is,  we  shall  determine  which  types  of  gravity  anomalies  may 
theoretically  be  used  in  making  geoidal  height  or  deflection  computations  both 
in  the  classical  and  in  the  new  geodetic  theory.  We  shall  not  go  into  great 
detail  on  the  mechanics  of  computation  of  anomalies  since  this  has  been  covered 
in  many  places  but  shall  be  primarily  concerned  with  the  logic  related  in  the  use 
of  the  different  types  of  anomalies. 

CLASSICAL  THEORY 

From  the  classical  geodetic  viewpoint  the  problem  is  to  use  gravity 
observed  on  the  physical  surface  of  the  earth  which  is  not  an  equipotential  sur¬ 
face  in  conjunction  with  a  theory  which  relates  gravity  on  a  bounding  equipoten¬ 
tial  surface  with  its  shape  and  the  direction  of  the  force  of  gravity  on  it.  In 
this  case,  the  four  distinct  steps  mentioned  at  the  beginning  of  Section  4  are 
necessary.  We  shall  first  consider  the  types  of  anomalies  which  result  from 
carrying  out  the  first  two  steps. 

The  most  common  types  of  gravity  reduction  methods  are: 

1.  Simple  free- air 

2.  Free-air  anomaly  with  condensation 

3.  Rudzki  inversion 

4.  Isostatic  based  on  hypothesis  of 

a.  Pratt -Hayford 

b.  Airy-Heiskanen 

c.  Vening  Meinesz 

3 


We  shall  begin  with  the  isostatic  reduction  method  since  it  has  been 
so  often  used  and  is  a  clear  example  of  what  is  done.  Physically  any  type  of 
isostatic  anomaly  envisions  a  number  of  processes. 

From  the  value  of  gravity  measured  at  P  on  the  earth's  surface,  pro¬ 
ceed  as  follows: 

1.  Compute  the  effect,  on  gravity  at  P,  of  removing  all  mass  above  the 
geoid.  This  is  the  complete  correction  for  topography  around  the  world  extending 
above  sea  level. 

2.  Compute  the  effect,  at  P,  of  placing  this  mass  below  the  geoid  at 
some  point  say  between  30  and  50  km.  below  sea  level.  Ihe  above  corresponds  to 
the  reduction  for  topography  and  compensation  in  isostatic  theory.  However,  in 
the  isostatic  reduction  procedure,  mass  is  also  moved  from  below  the  ocean  so 

as  to  fill  the  ocean  with  matter  having  the  density  of  rock.  Note  that  this  mass 
transfer  is  not  required  in  order  to  use  the  theory  of  Section  3  since  the  oceans 
are  already  bounded  by  the  desired  equipotential  surface.  As  long  as  the  proper 
computations  are  made  as  to  the  effect  of  such  an  internal  mass  transfer,  it  is 
of  course  permissible.  Ihe  same  may  be  said  for  geologic  corrections  below  sea 
level. 

If  a  normal  isostatic  correction  is  used,  the  effect  of  the  mass  transfer 
on  gravity  at  the  earth's  surface  is  given  by 

— B.C.  +  L.I.C,  +  D.T.A.C. 


where : 

B.C.  is  the  complete  Bouguer  correction  as  modified  by  Bullard  to 
extend  only  through  zone  0  and  is  considered  a  positive  number. 

L.I.C.  is  the  local  isostatic  reduction  out  through  zone  0. 

D.T.A.C.  is  the  combined  topographic -isostatic  correction  for  material 
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beyond  zone  0  as  this  zone  is  defined  by  the  U.S.  Coast  and  Geodetic  Survey 
and  may  be  either  positive  or  negative. 

Then  for  a  land  station  if  g  is  the  observed  gravity  at  a  point  P 
on  the  earth's  surface  and  g  is  the  modified  value  of  gravity  at  P,  we  have 
for  an  inland  station 

g  --  g  -  B.C.  *-  L.I.C.  +  D.T.A.C. 

where  for  the  station  the  effect  of  taking  away  the  above  sea  level  near  mass 
is  negative,  the  effect  of  putting  in  a  "compensating"  mass  below  sea  level  is 
positive  and  the  effect  of  the  distant  changes  may  be  either  positive  or  negative 
depending  upon  the  location  of  the  station  with  respect  to  the  totality  of  the 
changes. 

3.  Compute  the  value  of  gravity  at  the  point  P'  on  the  surface  of  the 

—  dg 

geoid.  This  means  starting  with  g  at  the  point  P  and  using  along  the  normal 
between  P  and  P'  to  compute  the  value  of  gravity  at  P' .  The  actual  gradient  is 
not  known  so  the  average  ellipsoidal  gradient  is  normally  used,  i.e.,  we  apply 
the  "free -air"  correction  to  yield 

g'  =  g  +  F.A.C. 

4.  Note  the  effect  of  the  mass  transfer  on  the  equipotential  surface 
which,  prior  to  the  mass  transfers,  coincided  with  the  geoid.  This  equipoten¬ 
tial  surface  has  been  moved  slightly  and  we  call  its  present  position  the  co-geoid. 

5.  Move  all  matter  lying  between  the  geoid  and  co-geoid  vertically  to 
the  surface  of  the  co-geoid.  We  assume  that  the  effect  on  gravity  at  P'  of  this 
move  is  negligible. 

6.  Use  the  gravity  gradient  to  compute  the  gravity  at  P"  on  the  co- 
geoid  which  we  shall  call  g’ '  from  the  value  of  gravity  at  P'  by  making  the  cor¬ 
rection  h'  (r^-)  where  h'  is  the  distance  between  the  geoid  and  co-geoid  and 

oh 


dh 
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is  the  gravity  gradient.  This  is  the  so-called  indirect  effect  so  we  have 

g"  =  g'  +  I.E. 

where  I.E.  =  indirect  effect  =  h'(|^) 

7.  We  can  now  get  the  gravity  anomalies  needed  for  Stokes  or 
Vening-Meinesz  Theorem  by  computing  Ag  =  g"  -  y^  where  y q  is  the  theoretical 
gravity  according  to  the  international  formula. 


One  should  note  at  this  point  that  Vening  Meinesz  at  step  5  rather 
than  transferring  the  mass  between  the  co-geoid  to  the  co-geoid  surface  assumes 
that  it  is  moved  to  infinity.  In  this  case  the  effect  of  the  mass  transfer  on 
gravity  is  not  negligible  but  must  be  computed.  The  remainder  of  the  steps  are 
as  outlined  above  in  6  and  J. 

The  details  of  carrying  out  isostatic  reductions  have  been  covered  in 
a  number  of  places  (See  for  example  Heiskanen  and  Vening  Meinesz,  1958)  and  will 
not  be  repeated  here.  Following  the  seven  steps  given  above,  the  formula  for 
any  type  of  isostatic  anomaly  can  be  written 

Agl  =  (Ob  -  B.C.  +  T.C.  +  R.C.  +  L.I.C.  +  D.T.C.  +  h  +  h'  )  -  Th 


where : 


B.C. 

T.C. 

R.C. 


L.I.C. 

D.T.C. 


^h' 

dh 


simple  Bouguer  slab  correction 
terrain  correction  through  zone  0 

Bullard  correction  to  remove  effect  of  material  beyond  zone  0 
and  curvature  correction 

correction  for  compensation  of  material  through  zone  0 

correction  for  topography  plus  compensation  beyond  zone  0 

free -air  correction  to  geoid  level  from  surface  level  using 
measured  heights 

free-air  correction  from  geoid  level  to  co-geoid  level 
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There  are  of  course  a  number  of  types  of  isostatic  anomalies.  From  the 
geodetic  point  of  view  these  simply  mean  that  the  mass  is  transferred  to  different 
points  inside  the  geoid  in  different  systems.  As  long  as  the  effect  of  these 
transfers  iB  accurately  computed  it  does  not  matter  much  which  one  iB  employed 
except  possibly,  as  will  be  discussed  later,  in  relation  to  prediction  of  the 
anomalies. 

Another  type  of  mass  transfer  which  is  permissible  and  can  be  used  is  to 
transfer  all  mass  above  geoid  level  to  the  geoid  surface  itself.  It  can  be  consid¬ 
ered  a  surface  density  on  the  geoid.  This  is  the  condensation  reduction.  It 
can  be  shown  that  such  a  mass  transfer  has  a  negligible  effect  on  the  shape  of 
the  geoid.  The  steps  in  using  this  method  are  then  reduced  to 

1.  Calculate  the  effect,  on  gravity  at  P,  of  moving  all  external 

mass  on  to  the  surface  of  the  geoid. 

2.  Use  h  i.e.,  a  free-air  reduction  to  compute  gravity  on 

the  geoid. 

The  development  of  this  theory  can  also  be  stated  in  the  following 
slightly  different  form:  Consider  a  point  on  the  earth’s  surface  P  and  a  cor¬ 
responding  point  on  the  geoid  P'  a  distance  h  below  P.  Assume  the  earth’s 
masses  above  sea  level  are  in  place.  Now  we  wish  to  arrive  at  a  formula  far  ||j- 
at  points  between  P  and  P’. 

We  begin  with  Poissons'  equation  for  a  rotating  body  of  density  p 


a2w  +  +  afw  .  .  4nkp  +  2u2 


a*  ay  a*  (5-i) 

If  we  have  the  z  axis  in  the  direction  of  the  outward  normal  to  the 


equipotential  surface 

a2W  _  _d_  (3W)  „  _  88 

a  z2  a  *  (a  * )  a  z 

(5-2) 


Th^n  equation  (5-l)  becomes 


||  -  +  i!»)  +  *„kp  .  2u2 

3xz  dy 


(5-3) 


But  it  can  be  shown  that 


2  2 
aw  +  d  w 


ax  ay 


/I  ,  1  X 
(-  +  -  )  6 
x  y 


(5-4) 


where  rY  and  r„  are  the  radii  of  curvature  of  the  equipotertial  surfaces  in  the 

x  y 

xz  and  yz  planes. 

Then  equation  (5-3)  becomes 


||  =  -  8  (i  +  i  )  +  4„kp  .  2J 


(5-5) 


Thus,  if  gp  is  the  value  of  gravity  at  the  point  P  on  the  earth’s  surface,  the 
value  gp'  at  the  point  P'  on  the  geoid  with  the  mass  above  sea  level  still 


in  place  is 


“  8p  -  I*1  (-  8  Of  +  i  )  +  4*kp  -  2wz]d-. 


(5-6) 


where  the  minus  in  front  of  the  integral  comes  from  fact  we  are  integrating 
in  negative  z  direction.  Thus 

V  =  8p  -  4nkM  +  J*1  l«  <i  +  i  >  -  2  <*i2]  dz 

p  x  y  (5-7) 

The  integral  on  the  R.H.S,  of  equation  (5-7)  is  a  free-air  correction 
in  the  strict  sense.  It  is  usual  to  assume  that  the  simple  free- air  correction  is 
sufficiently  accurate  and  write  g^,  ■  -  4rckM  +  .  09406h 


or,  if  we  write 


g  *  g  +  . 09406h 

o  p 

8p»  -  80  -  4*kM 


Thus  according  to  Green's  equivalent  layer  theorem  the  gravity  field 
outside  the  geoid  can  be  taken  as  equal  to  the  attraction  of  a  layer  of  density 
-  M  on  the  geoid  plus  the  attraction  of  material  outside  the  geoid.  We 
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will  now  remove  the  material  outside  the  geoid  and  the  surface  density  of  -M 
and  compute  their  effects  on  gp,.  The  effects  on  the  shape  of  the  equipotential 
surface  passing  through  P'  will  be  negligible.  We  now  have  a  gravity  value, 
gp,  that  can  be  used  to  obtain  an  anomaly  Ag  =  g^t  "  Yq  which  may  be  substituted 
into  Stokes'  or  Vening  Meinesz'  formulae.  If  the  effect  of  removing  the  M  layer 
on  the  geoid  and  the  material  above  the  geoid  is  the  same  and  each  equal  to  +2nkM 
then  gp»  =  gQ  =  +  .09^06h  and  the  simple  free-air  anomaly  could  be  used.  The 

assumption  that  the  removal  of  the  -M  layer  and  the  material  above  sea  level  is 
equal  to  4TTkM  is  fairly  accurate  in  level  areas  since  it  simply  amounts  to  the 
assumption  that  both  can  be  treated  as  infinite  horizontal  sheets. 

It  is  not  very  accurate  in  areas  of  considerable  relief  since  the 
effect,  at  F',  of  condensing  the  actual  topography  onto  the  geoid  is  considerably 
different  from  4TTkM  and  a  condensation  reduction  must  be  applied.  We  can  return 
to  this  point  later  in  connection  with  use  of  the  free  air  anomaly  in  the  modern 
theory  and  the  reduction  of  deflections  of  the  vertical  to  sea  level  or  the  sea 
level  deflections  to  surface  level. 

From  the  point  of  view  of  this  section  the  Model  Earth  anomaly  is  an 
offshoot  of  the  condensation  reduction.  One  first  carries  out  a  lateral  transfer 
of  material  at  the  earth's  surface  so  as  to  smooth  out  the  earth's  topography. 

The  assumption  that  the  M  layer  on  the  geoid  and  the  material  above  the  geoid 
have  an  effect  4nkM  when  carrying  out  a  condensation  for  the  smoothed  earth  can 
then  be  made  everywhere  with  sufficient  accuracy.  Since  the  primary  purpose  of 
the  Model  Earth  method  is  to  produce  an  anomaly  which  is  slowly  varying  for 
interpolation  and  extrapolation  purposes  we  shall  defer  additional  discussion 
of  it  until  Section  6.  For  present  purposes  it  sinrply  amounts  to  another  way  of 
transferring  material  from  above  to  below  the  geoid. 
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The  other  type  of  gravity  anomalies  sometimes  used  are  the  inversion 
anomalies  commonly  known  as  the  Rudzki  anomalies  after  the  Polish  scientist  who 
proposed  the  method.  The  reason  for  this  reduction  was  the  desire  to  avoid  any 
deformation  whatsoever  of  the  sea  level  equipotential  (geoid)  while  carrying  out 
the  mass  transfers.  The  Rudzki  method  has  been  little  used.  The  equations  for 
its  computation  can  be  found  in  Tengstrom  (1962).  Strictly  from  the  point  of 
view  of  the  classical  theory,  the  Rudzki  reduction  leaves  something  to  be  desired 
as  it  does  not  actually  produce  a  model  where  all  matter  is  within  the  bounding 
equipotential  but  leaves  a  small  amount  of  matter  outside  this  surface. 

There  are  other  types  of  anomalies  which  could  be  used  in  the  classi¬ 
cal  theory  relating  to  different  kinds  of  mass  transfer  but  in  fact  none  have  been 
used  and  would  in  any  case  differ  only  in  computational  detail. 

The  choice  of  the  anomaly  which  is  best  to  use  in  the  classical 
theory  is  beset  by  many  difficulties.  Primarily  this  results  from  the  fact  that 
knowledge  of  density  distributions  above  sea  level,  vertical  gradients  above  sea 
level,  and  the  geometric  relation  between  points  on  the  earth's  surface  and  on 
the  geoid  are  assumed  known.  In  fact,  they  are  not  exactly  known  and  consider¬ 
ations  of  minimizing  errors  due  to  lack  of  knowledge  must  be  taken  into  account. 

The  important  points  to  be  remembered  from  the  above  discussion  are 
that  the  necessary  requisites  for  an  anomaly  for  use  in  the  classical  theory 
(considering  only  the  point  that  they  must  be  theoretically  capable  of  being  in¬ 
serted  in  the  classical  equations)  is  that  they  envision  transfer  of  all  mass  so 
that  it  lies  inside  a  bounding  equipotential  surface  and  then  computation  of  the 
force  of  gravity  on  this  surface,  using  observed  gravity  on  the  earth's  surface. 

To  do  this  requires  knowledge  of  densities  above  sea  level  to  insure 
proper  mass  transfer  and  of  densities  below  sea  level  to  obtain  the  proper  gravity 
gradient.  Strictly  from  this  viewpoint  one  needs  to  have  a  ’’geologic  correction” 
above  sea  level  but  not  below  sea  level.  Again  from  the  viewpoint  of  satisfying 
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the  conditions  for  using  the  equations  of  Section  3  it  should  be  pointed  out  that 
the  fact  that  iBOstatic  equilibrium  does  not  or  does  exist  is  entirely  irrelevent. 
It  is  important  only  from  the  point  of  view  of  interpolation  and  extrapolation 
of  anomalies. 

THE  NEW  THEORY 

Whereas  the  theoretical  development  of  the  Classical  Geodetic  Theory 
does  not  imply  the  use  of  any  particular  type  of  anomaly,  the  New  Theory  (see 
Section  4)  naturally  leads  to  the  free-air  anomaly.  The  question  arises  as  to 
whether  or  not  some  other  type  of  anomaly  might  be  used.  The  following  develop¬ 
ment,  based  on  the  works  of  Arnold  (19^1),  shows  that  all  of  the  types  of 
anomalies  used  in  the  classical  theory  can  be  used  in  the  new  theory.  It  also 
shows  that,  considered  from  the  point  of  view  of  the  new  theory,  one  kind  of 
anomaly  gives  as  accurate  a  result  as  another  regardless  of  the  extent  of  know¬ 
ledge  of  internal  density  distributions.  This  leads  to  the  result  that  the  choice 
of  anomaly  is  dependent  on  ease  of  practical  computation  and  utility  in  extrapol¬ 
ation  and  interpolation  from  known  values. 

Perhaps  one  should  point  out  that  in  the  new  theory  the  equality  of  all 
types  of  anomalies  should  be  expected.  The  development  after  all  is  based  on 
small  deviations  from  a  known  model.  The  exact  form  of  the  model  could  not  be 
expected  to  control  the  accuracy  of  the  method. 

Consider  a  model  of  the  earth  consisting  of  the  international  ellipsoid 
model  plus  some  other  assumed  mass  distribution,  the  only  limitation  being  placed 
on  this  mass  distribution  being  that  it  must  lie  completely  within  the  physical 
surface  of  the  earth.  Let  the  potential  due  to  this  mass  distribution  be  desig¬ 
nated  by  A. 

Then 

W  =  potential  due  to  actual  earth 
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U  =  potential  due  to  international  ellipsoid  model 

A  =  potential  due  to  assumed  densities  in  theoretical  model  other 
than  those  of  international  ellipsoid  model. 

T'=W  -  (U  +  A )  =  difference  between  earth's  potential  and  model 
potential 

N  =  distance  from  P  to  Q 

Now  we  will  proceed  in  an  exactly  analogous  way  to  the  procedure  in 
Section  4.  Let  P  be  a  point  on  the  earth’s  surface  and  Q  the  corresponding 
point  on  the  equipotential  of  U  such  that  Uq  =  Wp. 

Now  by  definition 

V°PtAPfT’p  (5-8) 

Again  using  a  Taylor  expansion  of  U  about  Q  and  assuming  N  is  small 

enough  to  use  only  the  first  term  we  get 


(M) 

'■av’Q 


-  V 


Or  substituting  into  (5-8)  we  get 


(5-9) 


W  =  U  -  Y  N  +  A  +  T  1 P 
p  x  Q  P 


but  by  definition  Uq  =  Wp  so  we  get  from  (5-10) 


(5-10) 


H  =  --E-  +  AP_  (5-11) 

VQ  YQ 

This  equation  is  analogous  to  equation  4-12  of  Section  4.  The  term 
^  was,  in  the  older  theory  where  P  was  on  the  geoid  rather  than  the  earth’s 

YQ 

surface,  the  separation  between  the  geoid  and  the  co-geoid.  In  equation  (5-ll), 
since  we  supposedly  know  the  density  distribution  used  in  our  model,  we  can  com¬ 
pute  Ap  and  thus  in  (5-ll)  the  only  term  remaining  continues  to  be  the  potential 


T 


P* 
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Now  let  us  get  the  relation  between  gravity  anomalies  and  the  potential 
T'p.  From  the  definition  of  W#i.e.t  W  =  U  +  A  +  T  we  have 

,dw,  ,31),  .  ,3  A,  .  ,3T'. 


_  /O  U  \  .  /O  A  .  /0  f  v 

(T5>p  -  %n>p  +  <STi)p  +  %!>, 


P  (5-12) 

where  the  derivative  is  in  the  direction  of  the  normal  to  W  and  the  subscript  P 
indicates  that  the  derivatives  are  to  be  evaluated  at  the  point  P  .  Now 

au  au 

the  negative  of  the  observed  gravity.  Also  differs  from  ,  the  derivative 

on  a  V 

in  the  direction  of  the  normal  to  equal  U  surfaces,  only  by  a  quantity  which  is 
proportional  to  the  deflection  of  the  vertical  so  we  shall  assume  that  it  is 
accurate  enough  to  write 


(§J?) 

'31)% 


,3  IK 
(5v)p 


-  V. 


(5-14) 

Then  again  assuming  that  we  can  expand  y  around  the  point  Q  and  that  N  is  small 
enough  to  keep  only  the  first  term  we  get 

Y„ 


YQ  +  (I^Q  N  • 


(5-15) 


Taking  note  of  (5-13),  (5-1*0  and  (5-15)  we  get  for  (5-12) 


Or  rearranging  terms 


-  8. 


-  Vo  •  Oo»  V  o,  ♦  <?!>„ 


3  A , 


P  'Q 

Substituting  (5-H)  into  ( 5-l6)  gives 


(g  -  v  )  =  (^Y)  n  _  ,dAx  ,3T' 

d  Yo'  Wq  ^)p  -  (^) 


3  v  Q 

.3A 


(5-16) 


(g  -  v  )  -  (in  (HP  +  AP)  3T' 

VgD  Yo'  VdV;Q  1  v.  +  -  VTtOn  - 


3n'P  -  Wp 


(5-17) 


or  rearranging  terms 


[  (Sr 


p  -  V  +  Op  -  $>Q  $  -  <&>Q  ¥  -  <|f) 


<•3 Yi  Ap 


Q  Q  Yp 


'3"n'p 


(5-18) 
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The  expression  on  the  left  hand  side  of  equation  (5-18)  represents  some 

dA 

type  of  anomaly  such  as  an  isostatic  anomaly.  The  quantity  (-— )«  would  in  the 

on  y 

case  of  the  isostatic  anomaly  be  the  attraction  at  the  earth's  surface  of  topo¬ 
graphy  plus  compensation.  The  quantity  <Sv)q  %  is  analogous  to  the  indirect 

effect  in  classical  theory.  However,  the  term  (^3£)_^K.  must  be  evaluated  at  the 

^v'QYq 

earth's  surface  rather  than  at  geoid  level.  The  difference  is  probably  not  very 
great.  The  above  considerations  show  that  if  proper  computations  are  carried  out 
all  types  of  anomalies  will  yield  the  same  result. 

One  should  note  the  difference  in  the  viewpoint  adopted  in  arriving 
at  equation  (5-18)  which  holds  at  the  earth's  surface  and  that  adopted  in  arriv¬ 
ing  at  an  analogous  equation  in  classical  theory  which  holds  at  co-geoid  level. 

In  the  case  of  the  classical  theory  we  have  thought  in  terms  of  transfer  of 
actual  masses  of  the  earth  from  one  point  to  another,  which  implied  knowledge  of 
these  masses,  and  use  of  gravity  observed  at  one  level  to  compute  gravity  on 
another  level  which  implied  knowledge  of  actual  vertical  gradients.  In  the 
viewpoint  used  in  arriving  at  (5-18)  we  think  in  terms  of  the  density  distribu¬ 
tion  of  a  model  which  we  may  hope  approximates  the  actual  density  distribution 
of  the  earth  from  the  viewpoint  of  ease  of  interpolation.  However,  any  density 
deviations  between  the  actual  earth  and  the  model  do  not  introduce  errors  into 
the  result.  Also,  in  arriving  at  ( 5 -18 )  we  need  only  vertical  gradients  of  the 
gravity  field  of  the  models  which  can  be  computed  to  any  desired  degree  of 
accuracy  rather  than  the  unknown  vertical  gradients  of  the  actual  earth.  Thus 
even  if  we  are  using,  for  example,  an  isostatic  anomaly,  the  new  theory  pro¬ 
vides  a  much  sounder  basis  for  assessing  the  accuracy  to  be  expected  and  for 
determining  the  direction  to  be  taken  to  improve  accuracy. 

In  the  next  section  we  will  consider  the  questions  which  arise  when 
interpolation  and  extrapolation  of  anomalies  become  important. 
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It  is  in  connection  with  this  problem  of  anomaly  interpolation  that 
the  application  of  geology  and  geophysics  becomes  important. 


< 

I 
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SECTION  6 

Problems  of  Interpolation  of  Anomalies 

In  Sections  2  through  5  we  examined  the  mathematical  derivation  of 
the  formulae  relating  gravity  and  deflections  of  the  vertical.  Through  the 
developments  given  in  Section  5  we  see  that  in  the  "New  Geodetic  Theory”  the 
accuracy  with  which  one  may  obtain  deflections  is  entirely  independent  of  any 
knowledge  of  the  internal  density  distributions  of  the  earth  provided  there  is 
a  complete  knowledge  of  the  gravity  field.  Thus  one  can  never  hope  to  improve 
the  accuracy  of  the  theory  by  shifting  from  one  type  of  anomaly  to  another  or 
by  using  any  increased  knowledge  of  internal  density  distributions. 

The  choice  of  an  anomaly  to  use  in  deflection  computations  will  there¬ 
fore  be  controlled  by  two  considerations  which  arise  in  applying  the  theory  to 
actual  deflection  computations.  The  first  consideration  is  ease  of  computation. 
If  it  is  possible  to  decrease  the  computational  labor  by  choosing  one  type  of 
anomaly  over  another,  the  anomaly  requiring  the  least  labor  should,  of  course, 
be  chosen  provided  both  yield  answers  of  sufficient  accuracy.  The  second  con¬ 
sideration  is  accuracy  of  gravity  interpolation.  In  practice,  gravity  is  not 
known  everywhere  on  the  earth’s  surface  but  at  a  finite  number  of  points.  At 
intermediate  points  gravity  must  be  determined  by  some  form  of  interpolation. 

If  one  type  of  gravity  anomaly  can  be  more  accurately  interpolated  than  others 
it  should  improve  the  computed  deflections  by  providing  a  more  accurate  esti¬ 
mate  of  gravity  between  observation  points. 

In  arriving  at  a  final  choice  of  anomaly  an  attempt  must  be  made  to 
reconcile  the  two  considerations  so  that  the  anomaly  chosen  for  substitution  in 
the  deflection  formulae  yields  acceptable  accuracy  with  a  minimum  of  labor. 

This  means  that  one  must  determine  at  what  point  an  increase  in  accuracy  no 
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longer  justifies  the  increase  in  labor.  Another  point  to  be  kept  in  mind  in 
investigating  the  practical  utility  of  various  types  of  gravity  anomalies  is 
the  fact  that  one  type  of  anomaly  can  be  converted  to  another.  Thus  it  might 
be  possible  to  use  one  type  of  anomaly  for  interpolation  and  convert  to  another 
type  of  anomaly  for  substitution  into  the  deflection  formulae.  We  shall  first 
investigate  the  relative  accuracy  with  which  the  various  types  of  gravity  anom¬ 
alies  can  be  interpolated. 

Before  discussing  the  individual  types  of  anomalies  it  will  be  best 
to  examine  some  general  considerations  connected  with  interpolation  of  gravity 
anomalies.  The  ultimate  aim  in-so-far  as  interpolation  is  concerned  is  to  ob¬ 
tain  a  type  of  gravity  anomaly  which  varies  as  nearly  as  possible  linearly  be¬ 
tween  points  of  observation.  Using  this  type  of  anomaly  a  contour  map  prepared 
in  the  usual  way  would  have  maximum  accuracy.  An  alternative  would  be  to  use 
an  anomaly  which  might  deviate  from  linear  variation  between  observation  points 
but  in  such  a  way  that  the  deviations  could  be  predicted.  Such  a  procedure  might 
be  possible  for  example  by  use  of  geologic  and  geophysical  knowledge  to  predict 
non-linear  variations. 

In  order  for  a  gravity  anomaly  field  to  vary  nearly  linearly  between 
observation  points  it  must  contain  only  components  whose  wavelengths  are  much 
greater  than  the  distances  between  observation  points.  Any  gravity  anomaly  is 

difference  between  observed  gravity  and  the  gravity  field  of  a  density  model. 

Thus  to  produce  an  anomaly  field  which  contains  only  long  wavelength  components, 
the  density  model  used  in  anomaly  computation  must  correspond  as  nearly  as  possi¬ 
ble  to  the  actual  earth  with  respect  to  those  density  distributions  which  produce 
short  wavelength  gravity  components.  Long-wavelength  and  short-wavelength  are 
of  course  only  relative  as  used  here.  A  particular  wavelength  might  be  con¬ 
sidered  "long"  or  "short"  depending  upon  the  distance  between  observation  points. 
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It  is  clear  that  two  factors  control  the  width  of  the  gravity  effect  of  an 
anomalous  density  distribution.  These  are  the  horizontal  dimensions  of  the 
anomalous  density  distribution  and  its  proximity  to  the  surface.  Thus  the  ' 
density  model  used  in  computing  a  gravity  anomaly  for  interpolation  must  corres¬ 
pond  to  the  density  of  the  actual  earth  in  so  far  as  density  variations  of  less 
than  a  given  horizontal  dimension  and  given  depth  are  concerned.  The  magnitudes 
of  the  limiting  values  depend  upon  the  distance  between  observational  points. 

The  remarks  of  the  last  few  paragraphs  may  seem  to  be  simply  a  statement  of  the 
obvious  but  overlooking  this  obvious  can  lead  to  considerable  extraneous  dis¬ 
cussion. 

Let  us  now  turn  to  the  problem  at  hand — the  use  of  gravity  in  inter¬ 
polation  of  deflections  of  the  vertical.  The  accuracy  of  the  computed  deflec¬ 
tions  is  critically  dependent  upon  accurate  knowledge  of  the  gravity  field  in 
the  near  vicinity  of  the  deflection  station,  i.e.  within  100  kms.  Therefore, 
one  would  hardly  attempt  deflection  interpolation  without  a  fairly  dense  set  of 
observations.  As  a  maximum  for  the  distance  between  stations  30  km  is  a  reason¬ 
able  value.  Thus  the  problem  to  be  investigated  is:  What  density  distributions 
of  the  actual  earth  must  be  incorporated  into  a  density  model  in  order  to  pro¬ 
duce  an  anomaly  which  varies  smoothly  over  distances  of  30  km  or  less? 

By  far  the  most  important  near  surface  horizontal  density  variation 
producing  short  wavelength  gravity  effects  is  topography.  One  of  the  basic  re¬ 
quirements  for  any  smoothly  varying  gravity  anomaly  field  is  that  the  density 
model  used  in  its  computation  accurately  take  into  account  those  topographic 
masses  of  the  actual  earth  which  are  near  the  computation  stations. 

The  common  types  of  gravity  anomalies  which  result  in  a  smoothly  varying 
gravity  field  are: 

1. )  Bouguer  anomalies 

2. )  Isostatic  anomalies 

3. )  Model  earth  anomaly 

4. )  Rudzki  inversion  anomalies 
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In  all  of  these  types  of  anomalies  an  attempt  is  made  to  account  for  the  actual 
topography  of  the  earth  near  the  computation  point  in  the  density  model  used. 

In  the  following  discussion  we  shall  examine  the  ways  in  which  geologic  and 
geophysical  information  might  be  used  to  alter  the  density  models  so  as  to  pro¬ 
duce  an  even  more  smoothly  varying  gravity  anomaly  field. 

The  Bouguer  Anomaly 

If  the  Bouguer  Anomaly  is  to  be  used  for  interpolation  a  clear  under¬ 
standing  of  what  the  Bouguer  anomaly  implies  must  exist. 

The  normal  manner  of  writing  the  theoretical  model  attraction  for  the 
complete  Bouguer  anomaly  model  is: 

Theoretical  Gravity  of  Complete  Bouguer  Model  =  International 
Gravity  Formula  +  Attraction  of  Infinite  Slab  of  thickness 
equal  to  Height  of  station  -  Terrain  correction. 

In  abbreviated  form  this  becomes 

Th^  =  I.F .  2vkph  -  T.C. 

where  2irkph  is  the  well  known  infinite  slab  attratction. 

The  theoretical  model  attraction  for  the  normal  Bouguer  anomaly  models  seem 
irrational  in  terms  of  the  actual  earth.  A  closer  examination  shows  that  they 
are  not  quite  so  irrational  as  they  first  appear. 

First,  although  the  term  2xkph  theoretically  applies  to  an  infinite 
slab,  more  than  99$  of  this  attraction  results  from  the  part  of  the  slab  material 
lying  within  100  miles  (l66  km)  of  the  station.  Thus  to  a  close  approximation 
the  Complete  Bouguer  anomaly  represents  a  model  in  which  the  actual  earthte 
topography  out  to  1 66  km  has  been  superimposed  on  the  International  Ellipsoid 
model.  Bullard  has  suggested  a  slightly  modified  complete  Bouguer  model  to 
correct  for  the  minor  differences.  It  is  given  by 
ThcMB  =  I-F.+2irkph  -  T.C.  +  b 
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vhere  the  email  term  b  (only  3  to  4  mgls  at  most)  takes  care  of  the  difference 
between  the  inflate  slab  and  a  cylinder  of  radius  l66  km  and  also  corrects  for 
curvature  effects.  This  model  then  represents  exactly  the  International  ellip¬ 
soid  with  the  earth’s  topography  from  the  station  to  l66  km  superimposed. 

When  stated  as  above  it  becomes  apparent  that  the  density  model  used 


in  computing  a  Bouguer  anomaly  at  one  point  is  not  the  same  as  the  model  used 
at  another  point.  Perhaps  the  best  way  to  visualize  what  is  actually  indicated 
by  the  normal  Complete  Bouguer  Anomaly  is  as  follows.  Consider  a  density  model, 
M,  which  takes  into  account  the  actual  topography  above  sea  level  around  the 
world  and  the  mass  deficiencies  of  ocean  water.  We  might  write  the  attraction. 


where 


IF  =  attraction  of  International  ellipsoid  model 
NC  =  gravitational  effect  of  topography  within  166  km 
of  a  point 

DC  =  gravitational  effect  of  topography  beyond  l66  km 
of  a  point 


Both  of  the  functions  NC  and  DC  are  functions  of  position  on  the 
earth,  i.e.  of  the  coordinates  0,A  .  Thus  we  can  think  of  the  complete  Bouguer 
anomaly  an  being  written 

C.B.A.  =  gy  -  D.C.  (  e ,  X  ). 

The  function  D.C.  (  6  ,  X)  can  be  thought  of  here  as  simply  a  slowly 
varying  function  of  position  on  the  earth’s  surface.  The  Complete  Bouguer 
Anomaly  can  then  he  considered  as  simply  an  anomaly  computed  using  the  density 
model  M  plus  a  slowly  varying  function  of  position.  For  the  purposes  of  inter¬ 
polation  over  distances  of  30  km  or  less  the  function  D.C.  (  $  ,  A  )  will  be 
smoothly  varying  enough  to  cause  little  error  in  interpolation.  Of  much  greater 
importance  are  the  differences  between  the  model  density  and  the  actual  density 
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of  the  earth  within  l66  km  of  a  computation  point.  The  more  accurately  the 
model  densities  approach  the  actual  densities  of  the  earth  within  166  km  of  a 
point  the  more  smoothly  varying  will  be  the  resultant  anomaly.  The  alteration 
of  the  model  to  more  nearly  fit  the  actual  density  distribution  of  the  earth  as 
indicated  by  geologic  and  geophysical  information  will  be  termed  a  geologic 
correction. 

The  question  of  the  correct  density  assumption  often  arises  in  using 
the  Bouguer  anomaly.  The  answer  depends  upon  the  reason  for  computing  the 
Bouguer  anomaly. 

The  question  to  be  answered  in  the  present  investigation  is:  What 
density  should  be  used  for  the  Bouguer  correction  down  to  sea  level  and  for  the 
geologic  correction  below  sea  level  to  produce  the  most  easily  interpolated 
anomaly  when  observation  points  are  30  km  or  less  apart? 

Consider  first  the  suposition  that  one  computes  a  "geologically  cor¬ 
rected"  complete  Bouguer  anomaly  using  the  exact  densities  of  all  material 
above  the  depth  of  the  deepest  sedimentary  basin.  The  horizontal  changes  in 
this  anomaly  will  reflect  density  changes  at  depths  greater  than  the  depth  used 
in  computation  of  the  anomaly  plus  the  changes  in  the  effects  of  density  distri¬ 
butions  at  horizontal  distances  greater  than  l66  km.  Now,  if  the  surface  eleva¬ 
tion  is  5,000  ft.  and  the  bottom  of  the  deepest  basin  is  at  -25,000  ft.,  one 
will  be  at  least  30,000  ft.  above  any  remaining  density  variations.  This  dis¬ 
tance  factor  in  addition  to  the  fact  that  any  density  differences  at  this  depth 
will  normally  be  small  and  features  broad  would  be  expected  to  result  in  a  gravity 
effect  that  is  slowly  varying  in  a  horizontal  direction.  From  this  reasoning 
one  would  expect  this  type  of  anomaly  to  be  easily  interpolated. 

New  consider  an  anomaly  computed  using  a  density  of  2.88gm/cc  above 
sea  level  and  contrasting  density  below  sea  level  with  a  density  of  2.88gm/cc. 
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First  we  note  that  the  density  2.88gm/cc  arose  because  it  is  an  average  value 
for  the  mean  density  of  a  crustal  column.  In  the  first  part  of  section  IV  of 
the  final  report  of  contract  AF  23(601) -3*4-55  it  was  shown  that,  regardless  of 
the  density  of  the  part  of  the  crustal  column  above  sea  level,  if  one  assumes 
the  mean  crustal  density  remains  2.88gm/cc  and  mantle  density  remains  constant, 
one  would  get  the  correct  crustal  thickness  only  by  using  a  density  of  2.88gm/cc 
in  the  Bouguer  correction.  In  later  parts  of  section  IV  the  problem  of  changes 
in  mean  crustal  density  and  in  mantle  density  were  considered.  Let  us  tempo¬ 
rarily  assume  that  we  have  a  region  where  the  mean  crustal  density  remains 
2.88gm/cc  and  the  mantle  density  remains  constant.  The  question  is:  How  could 
the  reduction  of  the  observed  stations  in  this  area  using  a  density  of  2.88gm/cc 
in  the  Bouguer  anomaly  help  in  the  prediction  of  gravity  at  other  points. 

Let  us  consider  the  two  ways  in  which  the  2.88gm/cc  Bouguer  anomaly 
can  claim  to  be  better  than  the  Bouguer  anomaly  using  correct  densities  for 
interpolation  purposes .  These  ways  are 

1. )  The  2.88gm/cc  Bouguer  anomalies  themselves  could  be  more 

easily  interpolated. 

2. )  Seismic  information  could  be  used  to  predict  Bouguer 

anomalies  (  P  =  2.88gm/cc)  more  accurately  than  other 
types  of  Bouguer  anomalies. 

To  explore  the  first  contention  we  note  that  since  almost  all  topo¬ 
graphic  features  are  lower  in  density  than  2.88gm/cc  the  Bouguer  anomalies 
using  2.88gm/cc  will  retain  a  correlation  with  local  elevation  changes.  If  we 
rigorously  take  into  account  the  difference  in  density  between  material  present 
and  2.88gm/cc  either  above  or  below  sea  level  we  have  simply  added  a  constant 
to  the  result  we  would  have  gotten  using  any  other  density.  This  can  be  seen 
from  the  following  consideration. 

Let  V  be  the  volume  from  sea  level  to  a  depth  h  and  out  to  some  dis¬ 
tance  d  from  the  point.  Then  using  the  exact  densities  present  we  get  for  the 


attraction 


A  =  luj  r2  cos  v  dV 

When  using  the  density  difference  <y-  2.88  =  A<r  we  get  for  the  attraction 


difference 


A7  =  Is.J'AZ  cob<p  dV  =  k  J cos  *P  dV  -  k  2.88 


Thus  if  we  retain  the  same  area  of  integration  the  second  integral  is  a  constant 
at  all  points.  If  we  use  a  A <7  computation  only  for  sediments  and  not  for  the 
basement  we  will  introduce  spurious  effects  correlated  with  the  sediments.  It 
is  difficult  to  visualize  what  these  would  be  expected  to  represent. 

From  the  above  considerations  it  is  difficult  to  see  how  Bouguer 
anomalies  (  P  =  2.88)  would  be  more  easily  interpolated  than  those  using  the 
correct  densiJy  for  the  material  present.  Point  anomaly  values  computed  using 
P  =  2.88  would  include,  in  addition  to  the  effect  of  deep  seated  features, 
correlations  with  local  topography  and  near  surface  mass  distributions.  Examin¬ 
ation  of  the  two  types  of  Bouguer  anomalies  on  an  empirical  basis  also  failed  to 
show  any  interpolation  advantage  to  using  a  density  of  2.88gm/cc.  Figure  6-1 
shows  a  typical  profile  along  which  the  two  have  been  compared.  Additional  large 
scale  profiles  not  included  in  this  report  but  provided  to  ACIC  as  large  scale 
blackline  prints  show  similar  results. 


The  next  question  is:  Will  Bouguer  anomalies  P  =  2.88  averaged  over 
1°  x  1°  square  areas  by  more  closely  correlated  with  crustal  thickness?  If  this 
were  the  case  we  could  more  easily  estimate  these  anomalies  than  others  where 
there  were  seismic  measurements  of  the  earth’s  crust.  However,  this  is  not  a 
particular  reason  to  compute  Bouguer  anomalies  using  2.88gm/cc  in  an  area  where 
deflection  computations  are  to  be  undertaken  even  if  it  proves  to  be  true.  For 
deflection  interpolation  purposes  the  detail  needed  for  the  gravity  field  can 
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hardly  be  expected  to  be  provided  by  correlations  of  1  x  1  average  anomalies 
with  crustal  thickness. 

Perhaps  it  is  well  to  pause  here  and  again  emphasize  what  we  are  after. 
The  observed  gravity  at  the  earth’s  surface  is  a  rapidly  varying  function  pro¬ 
duced  by  all  of  the  masses  of  the  earth.  Because  points  where  gravity  is  ob¬ 
served  are  normally  too  far  apart  to  adequately  sample  the  short  period  changes 
accurate  interpolation  is  not  possible.  If  we  can  use  knowledge  of  existing 
density  distributions  to  compute  the  short  period  components  of  the  gravity 
field  and  remove  them  we  will  then  be  able  to  interpolate  the  long  period  por¬ 
tion  of  the  field  between  observation  points  with  adequate  accuracy.  If  in 
the  process  of  removing  the  short  period  components  of  the  gravity  field  we 
add  to  the  gravity  field  functions  which  have  only  long  period  components  the 
resultant  may  also  be  easily  interpolated. 

Two  factors  prevent  the  use  of  geologically  corrected  Bouguer  anomalies 
directly  in  deflection  computation  formulae.  First  as  noted  previously  the 
Bouguer  anomaly  as  usually  computed  does  not  use  the  same  density  model  at  all 
points.  To  utilize  Bouguer  anomalies  at  all  would  require  using  a  model  taking 
into  account  the  topography  around  the  world  with  the  increased  computations 
involved.  Second,  the  secondary  corrections  related  to  using  such  a  density 
model  would  be  long  and  complicated.  The  is  related  to  computation  of  the 
quantities  A  and  of  equation  (5-18)  for  this  model.  Also  one  often  has 
stations  which  are  sufficiently  close  together  to  allow  adequate  contouring  with¬ 
out  any  geologic  correction  in  one  area  while  in  an  adjacent  area  geology  must 
be  taken  into  account.  To  use  geologically  corrected  gravity  values  in  the 
formulae  would  require  that  the  effect  of  any  density  change  made  in  the  model 
be  computed  at  all  gravity  stations.  The  above  factors  combined  with  the  fact 
that  Pellinen  (1962)  has  developed  a  method  of  using  ordinary  complete  Bouguer 
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anomalies  directly  in  computations  leads  to  the  conclusion  that  the  best  method 
is  to  use  the  geologic  and  geophysical  information  for  interpolation  but  revert 
to  normal  complete  Bouguer  anomalies  to  draw  up  gravity  contour  maps  for  use  in 
the  actual  computations. 

Isostatic  Anomalies 

With  respect  to  the  problem  of  interpolation,  isostatic  anomalies  can 
be  thought  of  as  extensions  of  Bouguer  anomalies  just  as  the  Bouguer  anomalies 
are  themselves  an  extension  of  free-air  anomalies.  In  each  case  one  is  seeking 
to  improve  the  density  model  used  in  anomaly  computation  so  that  it  more  nearly 
approximates  the  actual  density  of  the  earth.  The  free -air  anomaly  is  computed 
using  as  a  density  model  an  ellipsoid  with  concentric  ellipsoidal  shells  of  con¬ 
stant  density.  In  computing  a  Bouguer  anomaly  this  density  model  is  altered  to 
take  into  account  the  earth’s  visable  topography  either  locally  as  is  usually 
done  or,  if  desired,  completely  around  the  world  using  a  constant  density.  In 
a  geologically  corrected  Bouguer  anomaly  lateral  variations  in  the  density  of 
the  upper  crust  as  indicated  by  geology  and  geophysics  are  included  in  the  den¬ 
sity  model  where  possible.  Both  gravity  observations  and  seismic  data  show 
that  large  scale  lateral  density  variations  occur  at  depth;  usually  in  such  a 
way  as  to  balance  the  topographic  and  other  near  surface  inequalities .  In 
isostatic  density  models  an  attempt  is  made  to  include  the  deep  seated  density 
variations  as  well  as  those  included  in  the  Bouguer  density  model. 

The  question  to  be  asked  here  is  whether  or  not  any  of  the  isostatic 
density  models  result  in  an  anomaly  which  is  more  easily  interpolated  over 
distances  of  30  km  or  less  than  is  the  geologically  corrected  Bouguer  anomaly. 
Both  from  the  theoretical  point  of  view  and  from  actual  observation  the  answer 
would  appear  to  be  no.  There  are  a  number  of  reasons  for  this.  First  note 
that  we  arp  only  interested  in  the  region  within  1 66  km  of  a  point.  That  the 
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effect  of  the  densities  assumed  in  the  isostatic  models  beyond  166  km  is  small 
is  obvious  from  the  publication  of  Karki,  et  al  (1961)  giving  the  effect  of  the 
distant  topography  plus  compensation  according  to  the  Airy-Heiskanen  model  T  =  30. 
Clearly  the  smoothly  varying  function  shown  there  could  be  of  no  effect  on  inter¬ 
polation  over  distances  of  30  km  or  less.  For  the  area  within  1 66  km  of  the  com¬ 
putation  point  the  amount  of  improvement  in  interpolation  accuracy  to  be  gained 
by  using  an  isostatic  anomaly  depends  upon  both  the  accuracy  with  which  the 
compensating  densities  of  the  model  actually  reflect  the  conditions  within  the 
earth,  and  the  degree  to  which  "local"  compensation  exists.  With  respect  to  the 
accuracy  of  the  density  model  assumed  it  is  becoming  increasingly  clear 
(Woollard,  1962)  that  compensation  is  not  achieved  in  a  uniform  manner  everywhere 
but  is  acnieved  by  complex  variations  in  crustal  thickness,  crustal  density  and 
upper  mantle  density.  Thus  a  density  model  which  assumes  that  compensation  is 
achieved  by  variation  of  a  single  parameter  such  as  crustal  thickness  is  not 
likely  to  give  more  than  a  rough  approximation  to  the  actual  density  distribu¬ 
tion.  The  fact  that  often  an  isostatic  density  model  will  nearly  explain  the 
observed  gravity  even  where  seismic  results  indicate  it  does  not  correspond  to 
the  actual  density  distribution  is  in  itself  an  indication  that  the  gravitational 
effect  of  the  compensating  masses  is  normally  a  rather  smoothly  varying  function. 

The  type  of  compensating  masses  most  likely  to  have  significant  short 
period  effects  in  the  gravity  field  would  be  local  compensation  for  features  of 
small  lateral  extent.  Every  indication  is,  however,  that  such  local  compensa¬ 
tion  does  not  exist.  The  finite  strength  of  the  earth* s  crust  would  tend  to 
make  local  compensation  unlikely  and  investigations  of  gravity  traverses  across 
features  of  small  lateral  extent  (Woollard,  1962)  indicate  that  local  compensa¬ 


tion  does  not  exist. 
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In  Figure  6-2  an  example  is  given  of  a  typical  profile  comparing 
isostatic  and  complete  Bouguer  anomaly.  It  can  be  6een  that  although  some 
large  scale,  long -wave length  trends  may  be  nearly  removed  by  using  an  isostatic 
density  model,  variations  over  distances  of  significance  for  interpolation  30  km 
or  less  are  virtually  unchanged.  The  same  result  may  be  seen  on  the  large  scale 
profiles  were  also  used  in  the  investigations  and  provided  to  ACIC  as  blackline 
prints.  The  general  conclusion  is  that  for  interpolation  over  short  distances 
of  30  km  or  less  not  enough  additional  interpolation  accuracy  would  be  attained 
to  justify  the  extra  work  entailed  in  isostatic  anomaly  computation  and  deter¬ 
mination  of  the  secondary  corrections  where  the  isostatic  anomalies  were  used 
directly  in  deflection  computations. 

Model  Earth  Anomalies 

In  so  far  as  interpolation  is  concerned  the  utility  of  the  Model  Earth 
Anomaly  can  be  quickly  determined  from  the  formula  for  its  computation.  The 

Model  Earth  Anomaly  can  be  written 

agme  =  g-  t"l'-T-A(H-h)]  =  g  -  (7+Ah-T)+AH 

•where  &%me=  Model  Earth  Anomaly 
g  =  observed  gravity 

7  =  International  formula  value  plus  free-air  correction 
T  =  terrain  correction 
A  =  Bouguer  correction  factor 
h  =  actual  height  of  station 
H  =  height  of  Model  Earth  at  station 

From  the  meaning  of  the  symbols  we  see  that  the  Model  Earth  Anomaly 
can  be  written 

Ag  =  Complete  Bouguer  Anomaly  +-  AH  (6-2^ 

me 

Since  the  height  H  is  the  weighted  mean  height  of  an  area  of  100  km 
or  more  in  radius  around  a  point  it  must  be  a  very  slowly  varying  function  of 
position.  Thus  the  model  earth  anomaly  will  be  no  better  than  the  Complete 
Bouguer  Anomaly  for  interpolation  over  short  distances  of  30  km  or  less  such  as 
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ve  are  interested  in  here.  Since  the  Model  Earth  Anomaly  will  be  no  better 
for  interpolation  than  the  complete  Bouguer  anomaly  and  it  requires  more  compu¬ 
tation  to  attain  this  anomaly  it  does  not  appear  to  have  any  particular  value 
for  deflection  interpolation  computations. 

Before  leaving  the  Model  Earth  Anomaly  it  is  perhaps  worth  while  to 
consider  the  reason  for  its  developement  since  it  so  well  illustrates  the  dif¬ 
ferences  in  outlook  between  the  Classical  and  the  New  Geodetic  Theory. 

The  Model  Earth  Anomaly  was  developed  for  use  with  the  classical 
theory.  Its  aim  was  to  develop  an  anomaly  which  could  be  interpolated  and 
extrapolated  as  easily  as  the  isostatic  anomaly,  would  be  more  easily  computed 
than  the  isostatic  anomaly,  would  avoid  assumptions  concerning  densities  down 
to  sea  level,  and  would  produce  negligable  deformation  of  the  geoid. 

With  the  advent  of  the  ’’New  Theory”  the  questions  of  incorrect  density 
assumptions  and  deformations  of  the  geoid  are  no  longer  important.  Thus  the 
Model  Earth  Anomaly  has  lost  much  of  its  reason  for  existing.  The  Model  Earth 
Anomaly  is  not  appreciably  easier  to  compute  than  an  Isostatic  Anomaly  in 
cases  where  tables  such  as  those  of  Karki,  et  al  (1961)  are  available  for  de¬ 
termining  the  distant  topographic -isostatic  effect.  In  this  case  the  primary 
computational  labor  in  either  method  is  related  to  making  the  terrain  correction 
out  to  166  km. 

Since,  as  we  have  seen,  the  Model  Earth  Anomaly  is  simply  the  complete 
Bouguer  anomaly  with  a  mean  elevation  term  added  it  should  be  smoothly  varying. 
However,  DeGraff  Hunter^  contention  that  the  Model  Earth  Anomaly  is  normally 
somewhat  smaller  in  magnitude  and  smoother  than  the  isostatic  anomaly  seems 
theoretically  somewhat  unlikely.  Since  the  Model  Earth  Anomaly  is  very  similar 
to  an  average  free-air  anomaly  and  average  free-air  anomalies  are  positively 
correlated  with  average  elevation,  we  would  theoretically  expect  the  Model  Earth 
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Anomaly  to  run  consistently  positive  in  mountain  areas.  Figure  6-3  illustrates 
this  tendency  in  an  area  in  the  Rocky  Mountains. 

Recently  Saunders  (1963)  has  made  a  study  of  Model  Earth  Anomalies 
in  which  he  computed  these  anomalies  for  a  large  area  in  Colorado.  His  re¬ 
sults  also  show  the  Model  Earth  Anomaly  to  be  consistently  more  positive  than 
the  isostatic  anomaly  in  mountainous  areas. 

Rudzki  Anomalies 

The  relative  accuracy  of  Rudzki  anomalies  for  interpolation  can  be 
seen  from  the  formula  for  the  attraction  of  the  Rudzki  density  model.  This 
formula  can  be  written 

g^  sb  I.F.  +  Top  -  T.M. 

where 

I.F.  =  International  formula  attraction 

Top  =  Attraction  of  topography  around  the  world 

T.M.  =  Attarction  of  deviation  masses  below  sea  level 

This  attraction  formula  resembles  that  of  an  isostatic  anomaly . 

However,  in  the  case  of  most  isostatic  models  the  deviation  mass  is  chosen  so 
as  to  try  to  approximate  the  actual  density  distributions  below  sea  level.  In 
the  Rudzki  density  model  the  deviation  mass  is  chosen  in  order  to  avoid  secondary 
corrections  when  using  the  classical  theory.  Clearly  an  anomaly  produced  by 
using  a  density  model  in  which  a  portion  of  the  density  distribution  does  not 
attempt  to  correspond  to  that  of  the  actual  earth  should  not  be  expected  to  be 
expecially  useful  for  interpolation  in  general.  The  Rudzki  anomalies  differ 
from  the  Bouguer  and  isostatic  anomalies  only  with  respect  to  the  long  period 
components  of  the  anomaly  field.  Thus  it  would  be  possible  to  apply  geologic 
corrections  to  a  Rudzki  anomaly  and  have  a  suitable  anomaly  for  inte rpolation 
over  short  distances.  Since,  however,  they  would  be  no  better  for  interpolation 
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than  Bouguer  anomalies  they  have  no  theoretical  advantage  in  the  "new  theory” 
and  would  require  much  greater  amounts  of  labor  to  compute.  They  have  not  been 
cons idered. 

Summary 

At  the  risk  of  being  repetitious,  the  results  of  this  section  will 
be  summarized  here  since  the  approach  used  and  the  conclusions  reached  vary 
somewhat  from  those  anticipated  when  the  investigation  wan  undertaken.  In  the 
contract  specifications,  profiles  from  ”Geoid  Correlations:  Test  Profiles” 
were  designated  as  profiles  to  be  studied  in  order  to  '’determine  the  best  geo¬ 
physical  reduction  method  to  transform  observed  gravity  into  anomalies  of  geo¬ 
detic  accuracy  and  significance".  At  the  time  this  study  was  begun  it  was  clear 
that  the  use  of  geologic  and  geophysical  information  could  improve  the  accuracy 
of  interpolation  and  extrapolation  of  gravity  values.  However,  it  was  not  clear 
in  what  way  this  fact  could  be  incorporated  into  the  theory  of  physical  geodesy 
in  order  to  achieve  the  most  improvement  in  the  computed  geodetic  quantities. 
Primarily  this  resulted  from  the  confusing  claims  and  counter  claims  of  propo¬ 
nents  of  the  various  types  of  anomalies.  Thus  to  begin  this  study  geodetic 
theory  was  carefully  examined.  As  a  result  of  this  part  of  the  study  it  became 
clear  that  modern  geodetic  theory  did  not  indicate  any  theoretical  advantage 
for  one  type  of  anomaly  over  another.  Most  of  the  apparent  advantages  and  the 
theoretical  problems  which  lead  to  them  no  longer  existed  when  using  the  "New 
Geodetic  Theory".  Instead,  the  accuracy  of  the  results  obtained  with  the  "new” 
theory  were  dependent  primarily  on  the  accuracy  with  which  the  gravity  field 
could  be  interpolated  between  observation  points.  The  result  of  this  is  a  re¬ 
orientation  in  the  way  gravity  anomalies  are  viewed.  Previously,  one  considered 
an  anomaly  to  result  from  computing  the  effects  on  the  observed  gravity  of 
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transfer  of  mass,  and  moving  from  one  point  to  another  in  space  before  comparison 
with  the  observed  gravity.  In  the  present  situation  one  can  think  of  the  gravity 
anomaly  as  simply  the  result  of  comparing  observed  gravity  at  the  point  of  obser¬ 
vation  with  the  gravity  predicted  at  that  point  by  some  density  model.  Thus  an 
anomaly  is  simply  the  difference  between  the  actual  density  distribution  of  the 
earth  and  the  density  distribution  of  the  model.  Viewed  in  this  way,  one  can 
almost  decide  even  before  studying  any  profiles  what  type  of  anomaly  will  be  best 
for  interpolation.  The  more  nearly  a  density  model  approximates  the  density 
distributions  of  the  actual  earth,  the  more  accurately  the  anomaly  computed  from 
it  can  be  interpolated. 

For  the  particular  application  under  consideration  in  connection  with 
this  contract  —  interpolation  of  deflections  of  the  vertical  --  a  reasonable 
amount  of  gravity  data  is  absolutely  necessary  to  achieve  an  answer  of  satisfac¬ 
tory  accuracy.  The  question  of  interpolation  over  short  distances  was  therefore 
the  question  of  greatest  importance  and  the  one  most  closely  studied.  The 
anomaly  profiles  prepared  as  required  under  sections  3b  and  3c  of  the 
technical  specifications  of  the  contract,  bore  out  the  general  conclusions 
reached  theoretically  concerning  the  general  problem  of  interpolation  but  were 
too  gross  in  scale  to  be  of  great  use  in  the  question  of  interpolation  over  short 
distances.  These  profiles  are  being  provided  to  ACIC  at  original  scale.  They 
are  not  included  in  the  report  since  it  was  felt  that  after  the  reduction  nec¬ 
essary  to  obtain  a  reasonable  size  for  binding  with  the  report  they  would  provide 
no  useful  information. 

In  order  to  study  empirically  the  relative  accuracy  of  anomalies  for 
interpolation  over  short  distances,  a  number  of  comparative  profiles  were  pre¬ 
pared  using  closer  spaced  observational  data.  These  bore  out  the  theoretical 


expectations . 


-79- 


Of  all  the  type6  of  anomalies  computed  without  using  geologic  and 
geophysical  control,  the  complete  Bouguer  anomaly,  using  as  nearly  as  possible 
the  actual  density  of  the  topography,  is  as  good  as  any  other  and  requires  less 
computation.  Anomalies  such  as  isostatic  and  Rudzki  anomalies  differ  from  the 
complete  Bouguer  anomaly  only  in  ways  that  do  not  effect  interpolation  over 
short  distances.  No  special  method  of  utilizing  geologic  and  geophysical  data 
was  discovered  under  this  contract.  Actually,  none  could  be  expected.  Once  it 
was  established  that  the  problem  was  one  of  simple  interpolation  it  becomes  a 
question  of  applying  all  of  the  results  from  Woollard  (1962)  and  the  results  thuB 
far  determined  under  contract  A P  23  (60l)  -  3^79  with  ACIC  to  any  particular 
problem  at  hand.  If  actual  subsurface  densities  are  known,  an  altered  density 
model  to  take  these  into  account  may  be  constructed  and  anomalies  computed  using 
this  model.  Alternately,  empirical  correlations  may  be  used  between  geologic  and 
geophysical  quantities  and  various  types  of  anomalies. 

The  existence  of  empirical  methods  which  can  be  applied  to  improve 
interpolation  even  when  no  exact  density  model  can  be  determined  suggest  the  use 
of  geologic  and  geophysical  knowledge  to  improve  the  interpolations  of  some  type 
of  existing  anomaly  rather  than  to  compute  an  anomaly  in  which  the  geologic  and 
geophysical  information  has  been  inserted  directly  into  the  density  models. 

The  complete  Bouguer  anomaly,  with  geology  and  geophysics  utilized  for  interpola¬ 
tion  between  observation  points,  appears  to  be  the  best  type  of  anomaly  to  use  to 
produce  a  contour  map  for  use  in  deflection  interpolation  computations.  The 
complete  Bouguer  anomaly  also  has  advantages  over  other  types  of  anomalies  insofar 
as  ease  of  computation  of  the  anomalies  is  concerned.  In  the  next  section  the 
manner  in  which  the  theory  of  Sections  4  and  5  may  be  combined  to  allow  the 
Complete  Bouguer  Anomalies  to  be  easily  utilized  for  computation  is  presented. 


-80- 


Section  7 

Formulation  for  Solution 


The  problem  of  the  formal  setup  for  solution  of  the  deflection  compu¬ 
tations  is  important  from  the  point  of  view  of  practicality.  According  to  the 
formulation  in  the  new  theory  the  deflection  of  the  vertical  in  an  arbitrary 
direction  1  is 


e 


R  r 

4  tty  ^ 


(Ag  +  G  x  ) 


dS  d t 

If  dl 


dw 


Ag  2N.  dH 

(“  +  T>  dl 


(7-1) 


where 

S  =  Stokes  function 

♦  =  angle  between  radius  vectors  of  fixed  and  variable  point 
R^dw  =  incremental  surface  element  of  sphere  of  radius  R 
Y  =  theoretical  g  (international  Formula) 

A  or  =  free  air  anomaly 

ro 

H  =  normal  height  of  variable  point 
H0  a  normal  height  of  fixed  point 
h  =  H  -  Hq  . 

rQ=2  R  sin  -j  =  distance  between  fixed  and  variable  points 
N  =  height  anomaly 

Pellinen  (1962)  seeks  an  answer  by  making  use  of  the  formulation  given  in  section 
5  allowing  use  of  anomalies  other  than  free  air  in  the  equations.  Pellinen  uses 
the  complete  Bouguer  anomaly  in  the  sense  that  all  mass  above  sea  level  is  actu¬ 
ally  considered  around  the  world.  This  instead  of  a  Bouguer  Plate  corrected  for 
topography  is  theoretically  used.  The  following  development  is  taken  from  Pel¬ 
linen  (1962)  with  some  clarifying  remarks  added  and  misprints  corrected. 

To  arrive  at  the  formula  for  deflection  of  the  vertical  we  will  first 
consider  the  formula  for  height  anomaly.  Let  us  call  the  potential  due  to  the 
topography  T  and  the  potential  difference  between  the  potential  of  the  actual 
earth  and  the  model  (which  is  the  international  formula  plus  topography  above 
sea  level)  T*  . 
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Then 


W  -  U  +  T  +  T' 

P  P  P  p 


Then  according  to  section  5  ve  have 
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and 


where 


4  we  get 


A*'  ■  -  -  <§£> 

Ag-  -  [(8p  -  VQ)  +  (||)p 


(7-4) 


(7-5) 


Then  in  a  manner  completely  analgous  to  the  development  of  section 


T  *  ^  J*  g*  +  Gp  S  (cos  f)  dw 


where 


Then 


Gi  =  in  J 


R2  p  hAg' 


dw 


N' 
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(7-6) 

(7-7) 

(7-8) 


and 


N  -  Ib 


'q  (7-9) 

However,  there  can  he  no  zero  or  first  order  harmonics  in  N  therefore 
we  must  write 


- 


TP  •  ho  ■  TPl 


(7-10) 

where  T  and  T,  are  the  zero  and  first  order  harmonics  in  the  formal  expansion 
Po  pt 

of  the  potential  of  the  topography. 

Now  we  note  that  the  attraction  of  the  material  above  sea  level  at  a 


-82- 


point  P  is 


4nkCTH  -  kcrR  P  (i  _  i) 
o  «  '  r  r 


°  ro  r  (7-11) 

2  o  i /p 

where  r  =  (r  +  h  )  '  =  distance  from  fixed  to  moving  point  on  surface  taking 
into  account  topography.  We  note  in  passing  that  the  second  term  comes  from  the 
integration  of  the  attraction  of  a  number  of  line  elements  with  rQ  and  r  being  the 
distances  to  the  two  ends  of  the  line.  Thus  the  attraction  of  material  above  sea 
level  is  broken  down  into  the  attraction  of  a  spherical  shell  of  thickness  H0 
plus  deviations  of  the  actual  topography  around  the  world  from  the  elevation  HQ. 
Now  going  back  to  equation  (7-5)  we  write 


[*•  ♦  n  ♦  «H  *  A>] 


But  according  to  (7-ll)  with  an  accuracy  of  the  order  of  Jj- 


4itkoH  -  kC7R2  P  (i  _  I)  dw 
o  °  r  r 


(7-12) 


(7-13) 


Ti  «  2*kCT«o  +  V2 


where  =  potential  of  second  term  on  RHS  of  (7-13) 

*n  neglecting  which  we  shall  investigate  later. 


Ag'  =  fig  +  \  |  -  2nkoH  +  Ag  ] 


(7-1*0 


(7-15) 


where  we  are  using  the  symbol  Ag^  to  represent  k  a  J  -  p-)dw.  Note  that 
Agp  is  essentially  a  terrain  correction  except  that  it  is  carried  completely 
around  the  earth.  Pellinen  states  that  it  actually  differs  very  little  from  the 
simple  terrain  correction.  Whether  or  not  we  accept  this,  the  difference  should 
contribute  little  when  used  for  interpolation  of  deflections. 

Now  we  shall  write  T  as  the  sum  of  two  parts^T®  the  potential  the  masses 
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of  the  earth  is  important.  Thus  we  might  write 


T„  -  kCTR2  r  dw  +  koR2  f  1  In  —  dw  +  V  ' 
P  J  r  *  o  r 

o  o 


(7-21) 


where  V*  represents  the  potential  of  topography  beyond  r^.  We  shall  assume  that 
this  can  be  ignored.  If  it  is  not  absolutely  negligable  in  magnitude  at  least 
it  is  slowly  varying  enough  to  be  ignored  in  deflection  interpolation  work.  Thus 


we  get 


AT  *  T  -  T° 
P  P 


-  koR2  J  dw  -  k  ct  R2  J  —  dw  +  kCTR2  J  dw 


2  p  h 


.2  p  H 


■jr 

+  kCTR2  Jo1  In  dw  OR  AT  .  kCTR2  J  (In  (r+h>  .  J>_) 


(7-22) 


or  expanding  ln(r  +  h)  in  series,  keeping  the  first  three  terms  and 

h  ~F°~ 

cancelling  the  terms  we  get 

P°  3  5 

AT  «  k  cr  R2  J  (-i  ---*  +  ™  --3 )  dw 


(7-23) 


Now  we  note  that  we  can  think  of  k  CT  H  in  (T-l6)  as  a  density  layer 

on  the  surface  of  a  sphere.  Then  introduce  a  function  0  0  associated  with  R  and 

P 

H  by  the  integral  equation 


G°  -  2 nkCTH  -  -2  kCTR  f  dw 
p  L  u  r 


(7-24) 


This  is  an  equation  exactly  analogous  to  equation  (4-51)  of  section  4 


and  the  solution  is 


k<j  J  dw  *  ZT51  -T  gp  8  (co8  dw 


(7-25) 


Or  from  equation  (7-16) 


-  K  -  -  £  x  gp  s(cos  ♦>  ** 


(7-26) 


where  the  zero  and  first  order  harmonics  have  been  left  out  in  the  solution. 
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Substituting  (7-l6)  into  (7-24) 


G°  =  2TTkaH  -  l 
P  * 


tS 


R 
—  o 


L-  J  (2«fcOH  -  |  2*)  dw 


(7-27) 

(7-23) 


Nov  going  back  to  equation  (7-3) 


N 


T  T°  AT 

2  -  P  4.  — 

VQ  ’  Yq  Yq 


(7-29) 


or  more  properly  since  there  can  be  no  zero  or  1st  order  terms  in  N 

At  -  At  -  At, 


T°  -  T°  -  T° 

P  PQ  pi 


or  substituting  from  (7-2 6) 

m  _  t>  0  o  m o  AT  -  Af  -  AT, 

j  (2itkcrH  -  |  *p)  S  (cos  *»r)  dw  +  - - ^-2 - i 

Q  v  yq 

Nov  substituting  (7-31)  and  (7-8)  into  (7-3)  ve  get 

R  P  -5  Tn  Af  -  AT 

N  =  7-—  J  (A- '  +  +  2jtkOH  -  \  S  (cos  iV)dw  +  c 

or  substituting  for  Ag'  from  (7-15)  and  collecting  terms 
K  “  zfy~  S  (AS  *  AZp  +  Sj  +  ^  -|2)  s  (cos  f,)  du  + 

Let  us  now  turn  our  attention  to  the  deflection  formulae  for 
type  of  development t  The  formula  for  deflection  is 

0  =  — —  (N )  =»  -g-  («~g  N  -f*  N) 


(7-30) 


(7-31) 


-  At 


yq 

(7- 

At  -  At 

o 

at 

O' 

>* 

! 

(7- 

T  T0 


T0 


this 


(7-34) 


or 


9  e  0  1  +  "0 


-1 

RYn 


(l-g  T '  (cp,  X,H)  +  Jg  T  (cp.X.H)} 


(7-35) 


but  as  shown  in  section  4  for  using  the  potential  where  it  is  known  we  have 

iv  T'  (9>x)  +  {^s'  +  (7-36) 

A  T  <t?«x’H>  "  II  <*-x>  *  {*s+T}|i  (7-37) 


but  substituting  from  (7-6) 


<7-  33) 


and  using  T  =  T°  +  AT  and  remembering  there  can  be  no  first  order  terms  we  get 


15  =  ol  {ye  <*°  -  t!>  +  h  <AT  -  Sf  i>  +  [*g  +  -R- 


3E) 


(7-39) 


but  from  (7-26) 


3  ■  £  f  •;  %  &  *> 


(7-40) 


Using  (7-38),  (7-39),  and  (7-40)  we  get 

0  =  e'  +  1  -  7I“-  |  <*3'  GX  +  G°p)  ||  ||  dw  +  A9p  - 


Ag  +T 


'1  311  _J 

J  3"S  Rvj 


(7-4X) 


where  A9p 


_ x  3  (At  -  ATt) 


V 


T0‘ 


or  since  Ag*  “  AS  +  |  |  -  2irkoH  +  Agp  and  G°  =  2jtl:OH  -  |  ~ 


9  =  tth  i*  <Ag  +  AsP  *  I  "I  +  Gl)  If  3T3  dw  *  A9p  -  [A=  +  x]  Il3  T^‘ 


Now  we  use  equation  (7-23)  to  "et 

-l-.oR  .a  ,  1  h3  .  3  h5 


A  8p 


vQ 


<- 


■5  — 3  ^  TO  — 5 

r„  r„ 

o  o 


)  dw 


(7-42) 

(7-43) 


But  one  can  write 


-37- 


Noting  that 


dr  ■  It  cos  ■»  d'yf 


dty  ■  cos  (rQ,  1)  d0 


(7-44) 


where  cos  (r0,  l)  is  the  cosine  of  the  angle  between  the  direction  in  which 
the  deflection  is  computed  and  the  direction  of  the  incremental  surface  ele¬ 
ment  from  the  station. 


Then  from  equations  (7-43)  and  (7-44)  we  get 


aeP  =  r  <■£  -|  cos  |  cos  (r0,i)  dw 

Yq  J  r  r 

but  dw  ■  sin  iJfd^dA  =*  2  sin  cos  -j\|fd^dA  =  ^  rdrdA 


,i  hJ  3  h5 


(7-45) 


h3  3  h 


J  <1  -3  -  o  cos  \  cos  drdA  (7-46) 


It  might  be  pointed  out  here  that,  instead  of  equation  (7-45)  one  has 
in  the  original  reference,  equation  19  of  that  reference  (Pellinen,  1962). 


4ep  B  kOTR  J  cos  ~  cos  (r  *  1)  (-  --  -  2  ll!)  dw 
2  0  2  r4  8  6 


(7-47) 


Clearly  the  factor  Yq  is  incorrectly  omitted  from  equation  (7~47).  A  question 
also  might  arise  with  respect  to  the  correct  power  of  R,  This  can  be  examined 
by  use  of  a  dimension  check.  The  units  of  k  are  cm^/gm  sec2.  Remembering  that 


dw  is  dimensionless  we  have  for  (7-45) 

on  o 

koR^  h  _  cm-3  gm 

” 7T  *  ~n+T  “  - 2  *  3  * 

‘Q  r  gmsec  cm 


cm 

n*T 


or  the  rhs  of  (7-45)  is  dimensionless  as  it  should  be.  If  we  use  equation  (7-47) 
with  Yq  added  to  the  denominator  the  rhs  would  be  which  would  be  in¬ 
correct  as  a  dimension. 

If  we  integrate  (7-46)  to  obtain  the  effect  of  an  area  lying  between 
A  a  A^,  and  A  =  A2  and  r  «  r-j_  and  r  =  rg  we  get  the  following,  using  the 


assumption  that  with  small  angles  cos 

A(A9p)  =  -12  cos  (rQ,l)  AA  {^-  <*g  -  -  |f-  (-^  -  )} 


rl  r2  rl  r2 


(7-40) 
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where  h  is  the  average  value  of  the  height  difference  and  A  A  =  -  Ag. 

For  practical  computations  the  above  equations  can  be  simplified  in 

the  following  manner.  The  quantities  3/2  and  Gi  beneath  the  integral 

R  1 

sign  in  equation  (7-42)  can  be  assumed  to  be  negligible  and  ignored.  Similarly 


the  term  I  Ag  + 


1  BH  i 

i  B6  Kvn 


can  be  assumed  to  be  negligible.  If  these 


assumptions  are  made,  equation  (7-42)  reduces  to 

6  "  J*  <As  +  As„>  I!  si  dw  +  a0p 


(7-49) 


Similarly  we  have  from  equation  (7-15) 


Ag1  =  Ag  +  Agp  -  2 jtIcaR 


(7-50) 


As  stated  previously  we  can  use  for  the  simple  terrain  correction.  In 

this  case  Ag'  is  the  complete  Bouguer  anomaly.  Thus  from  (7-50)  we  get 
Ag  4'  Agp  ■  Ag'  4-  2rtkaH  and  (7-49)  becomes 

e  =  7T^  I  (AS'  +  2*kCTH>  ||  ||  dw  +  4eP  (7-51) 

The  integral  on  the  rhs  of  (7-51)  is  exactly  analgous  to  the  normal 
integral  used  for  deflection  computation  in  classical  theory,  only  the  anomaly 
used  is  changed.  Thus  the  normal  Rice  template  may  be  used  to  solve  the  in¬ 
tegral  expression.  In  computations  it  is  found  best  to  use  the  Rice  template 
twice--once  for  each  integral  term. 

The  advantage  of  carrying  out  the  computations  in  this  way  is  that  the 
highly  variable  component  requires  estimation  of  average  elevations  rather  than 
average  anomalies  of  the  form  (Ag1  +  2 it  k  o  R).  Since  elevation  contour 
maps  are  already  available,  whereas,  the  special  anomaly  maps  required  would  be 
fully  as  complex  as  the  elevation  maps,  and  would  have  to  be  drafted  this  is  a 
tremendous  advantage. 


The  complete  Bouguer  anomaly  Ag'  on  the  other  hand  is  a  very  smoothly 
varying  function  and,  once  a  realistic  representation  of  the  gravity  field  has 
been  derived,  using  geological  and  geophysical  data,  maps  are  easily  prepared. 

An  additional  advantage  of  proceeding  as  in  the  last  paragraph  is 
that  the  term  A9p  requires  averaged  elevations.  Although  the  use  of  the  Rice 
template  for  integration  of  the  A9p  integral  does  not  produce  a  solution  in 
which  each  section  of  each  ring  has  the  same  weight,  the  average  elevations 
obtained  from  the  Rice  template  can  be  used  to  obtain  A0p  with  a  minimum  of 
additional  calculations.  Certainly  this  method  of  procedure  is  much  to  be  pre¬ 
ferred  over  going  through  the  entire  process  of  picking  average  elevations 
from  another  template. 

In  Part  2  of  the  Final  Report  on  this  investigation  the  procedure 
described  above  Trill  be  applied  to  actual  data  and  the  results  compared  with  ' 
observed  astro -geode tic  deflections.  A  detailed  outline  of  the  steps  followed 
in  applying  the  procedure  will  be  included. 
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